function R=heatimpl(N,M,K,T,ic);
% HEATIMPL  Solves heat eqn using backward in time
%    (implicit) finite differences.  Uses N by M (spatial 
%    by time) grid.  Uses Matlab's sparse and "A\b" to 
%    solve linear system. Returns stability ratio 
%    R=K dt/(dx)^2.  Specify initial condition in ic.
%    Boundary conditions fixed: u(t,0)=1, u(t,1)=0.
% Try "heatimpl(10,5,1,.2,zeros(1,9))" to solve problem 
%	  1.0.2.  Compare to "heatexpl(10,5,1,.2,zeros(1,9))".
% (Ed Bueler; 2/16/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(all(size(ic) == [1 N-1])), 
   error('init cond is wrong size'), end;
w(1,:)=[1 ic 0];

% create sparse tridiagonal matrix for implicit
e = ones(N-1,1);
B = spdiags([-R*e (1+2*R)*e -R*e],-1:1, N-1, N-1);
% use "figure(2); full(B)" to see B

for k=2:M+1
   w(k,1)=1; w(k,N+1)=0; % boundary conditions
   b=w(k-1,2:N)+[R zeros(1,N-2)]; % right side, = w_l + c
   w(k,2:N)=(B\b')'; % solves system
   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