clear;

mp = 50;  % number of terms for each infinite series

x = linspace(0,1,40);     % x-points in plot in x = 0..1
y = linspace(0,1,40);     % y-points in plot in x = 0..1
t = linspace(1/12,2,24);  % time points

% First compute e-values and coefficients
for n=1:mp
  for m=1:mp
    lambda(n,m) = (n*pi)^2 + (m*pi)^2;
    B(n,m) = (4/(pi*sqrt(lambda(n,m))))*((1/n)*(cos(3*pi*n/4) - cos(pi*n/4)) ...
	     *(1/m)*(cos(3*pi*m/4) - cos(pi*m/4)));
  end
end

% Compute full solution at all points in space and time
for k=1:length(t)
  for i=1:length(x)
    for j=1:length(y)
    
       u(i,j,k) = 0;
       
       for n=1:mp
	 for m=1:mp
	   u(i,j,k) = u(i,j,k) + B(n,m)*sin(n*pi*x(i))*sin(m*pi*y(j))* ...
		    sin(sqrt(lambda(n,m))*t(k));
	 end
       end
       
    end
  end
end

% Create the plots
mcount = 0;
for k=1:6
  figure(k)
  clf
 
  for m=1:4
    mcount = mcount + 1;
    subplot(2,2,m)
    mesh(x,y,u(:,:,mcount));
    colormap('jet');
    set(gca,'fontsize',10);
    title([ 'u(x,y,t) at time t = ',num2str(t(mcount))]);
    xlabel('x');
    ylabel('y');
    set(gca,'fontsize',10);
    axis([0 1 0 1 -0.75 0.75]);
  end
end

% Print the plots to JPG files
figure(1); print -djpeg95 memb1.jpg
figure(2); print -djpeg95 memb2.jpg
figure(3); print -djpeg95 memb3.jpg
figure(4); print -djpeg95 memb4.jpg
figure(5); print -djpeg95 memb5.jpg
figure(6); print -djpeg95 memb6.jpg