function U = crnich(f,g1,g2,a,b,c,n,m)
% U = crnich(f,g1,g2,a,b,c,n,m)
% Crank-Nicholson solution to the heat equation.
% This subroutine uses trisys.m to solve a tri-diagonal system.
% 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;
s1 = 2 + 2/r;
s2 = 2/r - 2;
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
Vd = s1*ones(1,n);
Vd(1) = 1;
Vd(n) = 1;
Va = -ones(1,n-1);
Va(n-1) = 0;
Vc = -ones(1,n-1);
Vc(1) = 0;
Vb(1) = feval(g1,k*0);
Vb(n) = feval(g2,k*0);
for j=2:m,
for i=2:(n-1),
Vb(i) = U(i-1,j-1) + U(i+1,j-1) + s2*U(i,j-1);
end
X = trisys(Va,Vd,Vc,Vb);
U(1:n,j) = X';
end