function lwfigure;
% LWFIGURE  Script (function) to reproduce figure 4.8 in Morton & Mayers 2nd ed.

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;   dt = dx;
   xs = x(2:J);  % x shortened
   xl = x(2:J) - dx/2;  xr = x(2:J) + dx/2;
   for k = 1:4
      U = exp(-10*(4*x-1).^2);
      for t = dt:dt:tf(k)
         U(1) = 0; % boundary condition
         U(J+1) = 0;  % extra boundary condition; L-W does not determine U(J+1)
         a = a_get(xs,t);
         al = a_get(xl,t);  ar = a_get(xr,t);
         at = at_get(xs,t);
         % use formula (4.44), Lax-Wendroff for linear, non-constant coeff;
         % L-W is FTCS + correction:
         dUdx = (U(3:J+1) - U(1:J-1)) / (2 * dx);
         diffus = ar .* (U(3:J+1) - U(2:J)) - al .* (U(2:J) - U(1:J-1));
         correction = (dt^2 / 2) * ( - at .* dUdx + a .* diffus / dx^2 );
         U(2:J) = U(2:J) - dt * a .* dUdx + correction;
      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

function z = a_get(x,t)  % accepts vector x and scalar t
z = (1 + x.^2) ./ (1 + 2*x*t + 2*x.^2 + x.^4);

function z = at_get(x,t)  % accepts vector x and scalar t
z = - ( 2*x .* (1 + x.^2) ) ./ ( (1 + 2*x*t + 2*x.^2 + x.^4).^2 );