% yn=coneq(y0,t0,tn,N)
% 
% Solve the convection equation on the periodic 
% domain with period 1 using the Lax-Wendroff
% method.  For more information see Morton and 
% Mayer Numerical Solution of Partial Differntial
% equations.
%
%    y0 -- initial distribution of heat
%    t0 -- time of the initial condition
%    tn -- ending time
%    N  -- number of time steps to take

function yn=coneq(y0,t0,tn,N)
    J=length(y0);
    dx=1/J;
    dt=(tn-t0)/N;
    a=1;
    nu=a*dt/dx;
    disp(sprintf('nu=%g',nu));
    yn=y0;
    for n=0:N-1
        ynp1(J)=1/2*nu*(1+nu)*yn(J-1)+...
            (1-nu^2)*yn(J)-1/2*nu*(1-nu)*yn(1);        
        ynp1(1)=1/2*nu*(1+nu)*yn(J)+...
            (1-nu^2)*yn(1)-1/2*nu*(1-nu)*yn(2);
        for j=2:J-1
            ynp1(j)=1/2*nu*(1+nu)*yn(j-1)+...
            (1-nu^2)*yn(j)-1/2*nu*(1-nu)*yn(j+1);
        end
        yn=ynp1;
    end
