clear;

%Number of wavelengths in the Box;
M=2;
%Nondimensional (maybe) advective speed C
C=1;
%Number of Periods
NT=10;
%Amplitude
A=1;

%Spatial step size and temporal step duration
dx=1;
dt=1;

%Calculate and show CFL
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).*dx;
Nmax=Xdomain.*NT./(C.*M);
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.*sign(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.
    else        %Evolve solution at later timesteps
        %Initialize temporary U(n+1) from U(n)
        U1L=UL;
        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);          
            %Unstable Euler Method
            %U1L(i)=UL(i)-(C.*UL(i+1)-C.*UL(i-1)).*dt./(2.*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 n+1 (Lax) 
        UL=U1L;
        UL(Imin+1)=U1L(Imax+1);
        UL(Imax+3)=U1L(2);
    end;
    %Show Results    
    plot(x,Ua,'k:',x,UL,'k-');
    axis([0 100 -2 2]);
    xlabel('Horizontal Position');
    ylabel('Amplitude');
    title(['Lax (Solid) and Exact (Dotted) Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    pause(0.005);
    
end;
