tnew=zeros(50,50)+10; 	% make a 100 x 100 plane, 10 degrees
told=tnew;		% as before, we need a told and tnew
time=0;
const=.1;
deltat=1.;


told(30:40,20:30)=80;  % make a small region 80 degrees

h=image(told,'CDataMapping','Direct')	% set up initial image
colormap('hot')		% pretty colors



while 1					% loop indefinitely
  time=time+deltat

% solve the equation!
% tnew is the next time step of the temperature, told is the last time
%    step, and tnew-told is the second derivative of the temperautre
%    (the OLD temperature) in the spatial direction

  for i=2:49 
    for j=2:49  % loop over both dimensions
      tnew(i,j)=told(i,j)+deltat*const*(told(i-1,j)+told(i+1,j)-2*told(i,j)+...
                                        told(i,j-1)+told(i,j+1)-2*told(i,j)) ;
                 
    end
  end

  told=tnew;	% make the old temp equal to the new temp to use during the
		% next go around

  set(h,'CData',tnew)		% display new temperature
  drawnow
end