function R=heatexpl(N,M,K,T,ic);
% HEATEXPL  Solves heat eqn using forward in time
%    (explicit) finite differences.  Uses N by M (spatial 
%    by time) grid.  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 "heatexpl(10,50,1,.2,zeros(1,9))" to solve problem 
%    1.0.2.  But note "heatexpl(10,40,1,.2,zeros(1,9))"
%    shows stability trouble and 
%    "heatexpl(10,30,1,.2,zeros(1,9))" blows up completely.
% (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];

for k=2:M+1
   w(k,1)=1; w(k,N+1)=0; % boundary conditions
   w(k,2:N)=w(k-1,2:N)+R*(w(k-1,1:N-1)-2*w(k-1,2:N)+w(k-1,3:N+1));
   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