function U = forwdif(f,g1,g2,a,b,c,n,m)
% U = forwdif(f,g1,g2,a,b,c,n,m)
% Finite difference solution to the heat equation.
% f is a function, input.
% g1 is a function, input.
% g2 is a function, input.
% a is the width of [0 a], input.
% b is the width of [0 b], input.
% c is the constant in the heat equation, input.
% n is the number of grid points over [0 a], input.
% m is the number of grid points over [0 b], input.
% U is the solution matrix, output.
h = a/(n-1);
k = b/(m-1);
r = c^2*k/h^2;
s = 1 - 2*r;
U = zeros(n,m);
for j = 1:m,
U(1,j) = feval(g1,k*(j-1));
U(n,j) = feval(g2,k*(j-1));
end
for i = 2:(n-1),
U(i,1) = feval(f,h*(i-1));
end
for j = 2:m,
for i = 2:(n-1),
U(i,j) = s*U(i,j-1) + r*(U(i-1,j-1) + U(i+1,j-1));
end
end