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