function U = dirich(f1,f2,f3,f4,a,b,h,tol,max1)
% U = dirich(f1,f2,f3,f4,a,b,h,tol,max1)
% Dirichlet solution to Laplace's equation.
% f1 is a boundary function, input.
% f2 is a boundary function, input.
% f3 is a boundary function, input.
% f4 is a boundary unction, input.
% a is the width of [0 a], input.
% b is the width of [0 b], input.
% h is the step size, input.
% tol is the tolerance, input.
% max1 is the maximum number of iterations, input.
% U is the solution matrix, output.
n = fix(a/h)+1;
m = fix(b/h)+1;
ave = (a*(feval(f1,0)+feval(f2,0)) ...
+ b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);
U = ave*ones(n,m);
for j=1:m,
U(1,j) = feval(f3,h*(j-1));
U(n,j) = feval(f4,h*(j-1));
end
for i=1:n,
U(i,1) = feval(f1,h*(i-1));
U(i,m) = feval(f2,h*(i-1));
end
U(1,1) = (U(1,2) + U(2,1))/2;
U(1,m) = (U(1,m-1) + U(2,m))/2;
U(n,1) = (U(n-1,1) + U(n,2))/2;
U(n,m) = (U(n-1,m) + U(n,m-1))/2;
w = 4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));
err = 1;
cnt = 0;
while ((err>tol)&(cnt<=max1))
err = 0;
for j=2:(m-1),
for i=2:(n-1),
relx = w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+ U(i-1,j)-4*U(i,j))/4;
U(i,j) = U(i,j) + relx;
if (err<=abs(relx)), err=abs(relx); end
end
end
cnt = cnt+1;
end