% Math 3150-1
% 
% Sample code for example 3.7.2:
%
% Heat equation on a square plate of side 1 with c=1/pi and
%    initial temperature distribution u(x,y,0) = f(x,y) = 100

% number of terms in the expansion
N=20; M=20;

% Time steps and final time
nt=20; Tmax=1;
ts = linspace(0,Tmax,nt+1);

% whether to take snapshots for several times
take_snapshots=1;
nsnaps   = 4;
snapbasename = '3_7_2';
snapcount = 0;
if (mod(nt,nsnaps)~=0) 
 error('number of snapshots must divide number of time steps');
end;


%%%% Code starts here

% setup a grid to plot the function
Nx=20; Ny=20; % number of points in x and y directions
x=linspace(0,1,Nx); y=linspace(0,1,Ny); 
[xx,yy]=meshgrid(x,y);

for it=1:length(ts),
 t = ts(it);
 ss = zeros(size(xx));
 for n=1:N,
  for m=1:M,
    % Sine Series coefficients for initial temperature distribution
    Amn = 400*(1-(-1)^m)*(1-(-1)^n)/pi^2/m/n;
    lambdamn = sqrt(m^2+n^2);
    ss = ss + sin(m*pi*xx).*sin(n*pi*yy)*Amn*exp(-lambdamn^2*t);
   end;
  end;
  % the rest of the loop handles plotting
  mesh(xx,yy,ss,'edgecolor','k','facecolor','none');
  axis([0 1 0 1 0 110]);
  title(sprintf('t = %d',t)); xlabel('x'); ylabel('y');
  pause(0.5); % 0.2 is the time we pause
  % take a snapshot if necessary
  if (mod(it-1,nt/nsnaps)==0 & take_snapshots)
   filename = sprintf('%s_%03d.png',snapbasename,snapcount);
   fprintf([ 'saving to file ' filename '\n']);
   print('-dpng','-r50',filename);
   snapcount=snapcount+1;
  end;
end;
