function xx=diffuse(period)
x=1:20; % x is the numbers 1 through 20
time=0;
const=.01;
tnew=zeros(1,20)+300; % initial conditions of 300 deg everywhere
told=tnew;% we are going to need two separate temperature vectors in
% order to simultaneously solve for the current temperature and
% not overwrite the previous temperature
h=plot(tnew); % a cool way to animate plots
xlabel('distance')
ylabel('temperature')
axis([1 20 250 350])
set(h,'EraseMode','none','LineWidth',3)
while 1 % a MATLAB command to loop until
% you press cntrl-C
time=time+1;
told(1)=300+30*sin(time/period*2*pi); % sinusoidally varying temperature
tnew(1)=told(1);
% at surface
% now 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:19
tnew(i)=told(i)+const*(told(i-1)+told(i+1)-2*told(i));
end
% that was it! those 3 lines were the full solution to the differential
% equation for a single time step.
told=tnew; % make the old temp equal to the new temp to use during the
% next go around
% plot(tnew,'r')
% plot(tnew,'b')
set(h,'YData',tnew,'Color',[0,0,1]) % overplot a blue line
set(h,'Color',[1,0,0]) % overplot a red line
drawnow % plot it!
end