function heat2d(n,delt,imp)
%HEAT2D  Demonstration finite differences for solution of 2D heat equation
%
% Demonstrate the forward/centered and backward/centered
% finite difference solution of the homogeneous heat equation
% with Dirichlet boundary conditions on the square [-1,1]x[-1,1].
%
% HEAT2D(n,delt[,imp])
%
% Uses an (n+1)x(n+1) grid including boundary points, so
% mesh size is h=2/n.  Timestep is delt.  Stability condition
% for forwarded/centered differences (the default) is
% delt/h^2 <= .25.  Compare n=20, delt=.002 (stable), with
% n=20, delta=.005 (unstable).
%
% If an optional third argument is given, then
% backward/centered differences are used (unconditionally stable).

if nargin > 2
  imp = 1;
elseif nargin < 2
  error('Wrong number of arguments')
else
  imp = 0;
end

% adjust n
n = round(max(n,3));
tmax = 1.;
h = 2/n;
[x,y] = meshgrid(-1:h:1);
G = numgrid('S',n+1);
A = -delsq(G);
sigma = delt/h^2;

% Initialize the solution
U = -.6*(3*(-x+y-1)/5 - (3*(-x+1.5*y)).^3 ).*exp(-(3*x).^2-(4*y).^2) ;
U(find(G==0)) = 0;

% u is a vector of the values of the solution at the interior mesh points
p = find(G>0);
u = U(p);


minu = min(U(:));
maxu = max(U(:));

set(gcf,'pos',[100 100 1.5*560 1.5*420])
colormap(jet)
surf(x,y,U)
caxis([-1 1])
axis off

axis([-1,1,-1,1,minu,maxu])
camzoom(1.75)
pause

for t = delt:delt:tmax
% advance one time step
  if imp == 0
    u = u + sigma*A*u;
  else
    u = (speye(size(A))-sigma*A)\u;
  end
  U = zeros(size(x));
  U(p) = u;
  surf(x,y,U)
  caxis([-1 1])
  axis off
  axis([-1,1,-1,1,minu,maxu])
  camzoom(1.75)
  pause
end
