% UPWINDFIGURE  Script to produce a figure similar to figure 4.8
% in Morton & Mayers, but using the upwind method.

JJ=[50 100];  % dx=[0.02 0.01];
tf=[0.0 0.1 0.5 1.0];
for s=1:2
   J=JJ(s);   dx=1/J;   x=0:dx:1;
   xs=x(2:J+1);  % x shortened
   dt=dx;   nu0=dt/dx;
   for k=1:4
      U=exp(-10*(4*x-1).^2);
      for t=dt:dt:tf(k)
         U(1)=0; % boundary condition
         a=(1+xs.^2)./(1+2*xs*t+2*xs.^2+xs.^4);
         nu=nu0*a;  %  nu is a vector (a function of x)
         U(2:J+1)=(1-nu).*U(2:J+1)+nu.*U(1:J);  % upwind method
      end
      Uex=exp(-10*(4*(x-tf(k)./(1+x.^2))-1).^2);  % exact soln
      subplot(4,2,2*(k-1)+s), plot(x,U,'.-',x,Uex,':')
      set(gca,'Visible','off')  % turn off axes
      if k==1,  text(0.2,1.2,['\Delta x = ' num2str(dx)]), end
      if s==1,  text(0.9,-0.3,['At t = ' num2str(tf(k))]), end
   end
end