function R=heatCN(N,M,K,T,ic);
% HEATCN: R=heatCN(N,M,K,T,ic)
%    Same calling sequence as heatexpl.m, heatimpl.m.
%    Solves heat equation using Crank--Nicholson finite 
%    diff method.  Returns stability ratio 
%    R=K dt/(dx)^2.  Uses tridiag.m to solve linear system.
%    (Easily rewritten to use "A\b".)
% Try "heatCN(10,5,1,.2,zeros(1,9))" to compare.
% Try "heatCN(100,50,1,.2,zeros(1,99))" to push it too far.
% (Ed Bueler; 2/27/02)

dx=1/N; dt=T/M; x=0:dx:1; t=0:dt:T;
R=(K*dt)/(dx^2); blowup=10.0;
w=zeros(M+1,N+1);
if not(isequal(size(ic),[1 N-1])) 
   error('init cond is wrong size'), end;
w(1,:)=[1 ic 0];

% create tridiagonal N-1 X N-1 matrix for left side
d = (1+R)*ones(1,N-1); % diagonal of A
a = -.5*R*[0 ones(1,N-2)]; % super diagonal
c = -.5*R*[ones(1,N-2) 0]; % sub diagonal

for k=2:M+1
   w(k,1)=1; % boundary condition
      for j=2:N % build right side
      b(j-1)=.5*R*w(k-1,j-1)+(1-R)*w(k-1,j)+.5*R*w(k-1,j+1);
   end;
   b(1)=b(1)+.5*R; % include boundary conditions
   b(N-1)=b(N-1)+0;
   w(k,2:N)=tridiag(a,d,c,b,N-1);
   % compare: w(k,2:N)=(B\b')';
   w(k,N+1)=0; % boundary condition
   if norm(w(k,:),inf)>blowup
      error(['blow up at step ' num2str(k)]), end;
end;

%display it
clf, mesh(x,t,w), %also try: contour(x,t,w)
xlabel 'x axis', ylabel 't axis', zlabel 'u axis',
view(37.5,30), axis([0 1 0 T 0 1.2]) % remove if using contour