clear;

%Number of wavelengths in the Box;
M=1;
%Nondimensional (maybe) advective speed C
C=0.9;  
%Number of Full-domain Periods to Run
NT=4;
%Wave Amplitude
A=1;

%Spatial step size and temporal step duration
dx=1;
%Uncomment to set dt for CFL=1
%dt=dx./C;
%Uncomment to set dt for 1, so that CFL=C (Comment for CFL=1)
dt=1;
%Calculate and show CFL number
CFL=C.*dt./dx

%Max domain size in space, and evolution duration time
Imin=0;
Imax=100;
Nmin=0;

%Max spatial extent and temporal duration
Xmin=dx.*(Imin-1);
Xmax=dx.*(Imax+1);
Xdomain=(Imax-Imin+1).*dx;
Nmax=Xdomain.*NT./C;
Tmin=dt.*Nmin;
Tmax=dt.*Nmax;


x=[Xmin:dx:Xmax]
%Open Figure
figure(1);

%Calculate later times t...
for n=Nmin:1:Nmax   
    %Define analytical solution, also used for initial value
    Ua=A.*sin(2*pi*(x./(Xdomain./M))-2*pi*((n.*dt)./(Xdomain./(C.*M))));
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        UL=Ua;  %Lax Init.
        ULW=Ua; %Lax-Wendroff Init.
    else        %Evolve solution at later timesteps
        %Initialize temporary U(n+1) from U(n)
        U1L=UL;
        U1LWa=ULW;
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1) for Lax Method
            U1L(i)=(1./2).*(UL(i+1)+UL(i-1))-(C.*UL(i+1)-C.*UL(i-1)).*dt./(2.*dx);
        end;
        for im=Imin-1:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1/2) for Lax-Wendroff 2-Step Method
            U1LWa(i)=(1./2).*(ULW(i+1)+ULW(i))-(C.*ULW(i+1)-C.*ULW(i)).*dt./(2.*dx);            
        end;
        %Perform second Lax-Wendroff step:
        %Update results for n+1/2 (Lax-Wendroff)
        ULWa=U1LWa; 
        %Initialize Array
        U1LW=ULWa;
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1) for Lax-Wendroff 2-Step Method
            U1LW(i)=ULW(i)-(C.*(ULWa(i)-ULWa(i-1))).*(dt./dx);            
        end;
        %Assign value at right boundary, such that domain is "periodic"
        %This means, physically, that whatever goes out on the right will 
        %come back in on the left.
        %Update boundaries for Lax 
        UL=U1L;
        UL(Imin+1)=U1L(Imax+2);
        UL(Imax+3)=U1L(2);
        %Update boundaries for Lax-Wendroff
        ULW=U1LW;
        ULW(Imin+1)=U1LW(Imax+2);
        ULW(Imax+3)=U1LW(2);      
    end;
    %Show Results    
    plot(x,Ua,'k:',x,ULW,'k--',x,UL,'k-');
    axis([0 100 -1 1]);
    xlabel('Horizontal Position');
    ylabel('Amplitude');
    title(['Lax (Solid), Lax-Wendroff(Dashed) and Exact (Dotted) Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    pause(0.005); 
    
end;
