% EXPLICITFIG  Script to make fig 2.2 in Morton&Mayer,2nd ed. ELB 2/4/07

J=20;  dx=1/J;  x=0:dx:1;  dt=[0.0012 0.0013];  tsteps=[0 1 25 50];
NN=406;  xx=0:.001:1;  % for exact soln use NN terms of Fourier sum
u0(1:J/2)=2*x(1:J/2);  u0(J/2+1:J+1)=2-2*x(J/2+1:J+1); % init cond; eqn(2.24)

for k=1:2  % k gives column in figure
   nu=dt(k)/(dx^2);  % each column has a nu value; see eqn (2.20)
   for l=1:4  % l gives row in figure
      % compute and plot approx solution
      U=u0;  % start over with initial condition
      for n=1:tsteps(l)  % does not execute if tsteps(l)<1
         % iterate equation (2.19):
         U(2:J)=U(2:J) + nu*( U(3:J+1) - 2*U(2:J) + U(1:J-1) ); end
      subplot(4,2,2*(l-1)+k),  plot(x,U,'.-')  % dots with lines
      % exact solution; put labeled axes on upper-left plot
      tt=tsteps(l)*dt(k);   uu=zeros(size(xx));
      sign=1;
      for m=1:2:NN  % see exercise 2.1 (page 56) and (2.11)
         c=m*pi; cc=c^2;
         uu=uu+sign*(8/cc)*exp(-cc*tt)*sin(c*xx);
         sign=-sign;  end
      hold on,  plot(xx,uu),  hold off
      if ((l==1)&(k==1))
         set(gca,'XTick',[]), set(gca,'YTick',[])  % remove ticks
         xlabel x, ylabel U
      else, axis off, end  % other plots simply get no axes
      if l==1, title(['\Delta t = ' num2str(dt(k))]), end
   end
end