clear;
%Boundary Conditions
%bc=0 for open, bc=1 for reflecting
bc=1;
%Number of waves in box (determines wavelength)
%For a good time, try a non-half-integer-multiple
M=4;
%Nondimensional (maybe) advective speed C
C=1;  
%Number of Wave Periods to Run
%And by periods, I mean full-domain propagation times!
NT=10;
%Wave Amplitude
A=1;
%Gaussian envelope "width" (stdev) as fraction of total runtime
G=1/6;
%Gaussian packet "start" time as fraction of total runtime
T=1/2;

%Spatial step size and temporal step duration
dx=1;
%Set CFL Number
%(Note that Lax will be quite diffusive for <1)
CFL=1;
dt=CFL.*dx./C;

%Max domain size in space, and evolution duration time
Imin=0;
Imax=100;
Xmin=-1; %One added for BCs
Xmax=101;
Xdomain=(Imax-Imin+1)*dx;
Nmin=0;
Nmax=Xdomain*NT/C;
Tmin=dt*Nmin;
Tmax=dt*Nmax;

%Domain X-axis Reference
x=[Xmin:dx:Xmax];

%Open Figure
figure(1);

%Calculate later times t...
for n=Nmin:1:Nmax   
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        %Define initial zero state
        Ua=x.*0;
        Va=x.*0;
        UL=Ua;  %Lax Init.
        VL=Va;
        ULW=Ua; %Lax-Wendroff Init.
        VLW=Va;       
        U1L=UL; %Initialize temporary U(n+1) 
        U1LWa=ULW;
        V1L=VL;
        V1LWa=VLW;
    else        %Evolve solution at later timesteps
               
        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*VL(i+1)-C*VL(i-1))*dt/(2*dx);
            V1L(i)=(1/2)*(VL(i+1)+VL(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*VLW(i+1)-C*VLW(i))*dt/(2*dx);       
            V1LWa(i)=(1/2)*(VLW(i+1)+VLW(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; 
        VLWa=V1LWa; 
        %Initialize Array
        U1LW=ULWa;
        V1LW=VLWa;
        
        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*(VLWa(i)-VLWa(i-1)))*(dt/dx);   
            V1LW(i)=VLW(i)-(C*(ULWa(i)-ULWa(i-1)))*(dt/dx);            
        end;
        
        %Set solution to updated values
        UL=U1L;
        VL=V1L;
        ULW=U1LW;
        VLW=V1LW;
       
        %Update boundaries for Lax 
        UL(Imin+1)=U1L(Imin+2)+A*sin(2*pi*((n*dt)/(Xdomain/(C*M))))*exp(-.5*((n*dt-Tmax*T)^2)/(Tmax*G)^2);
        UL(Imax+3)=U1L(Imax+2);
        VL(Imin+1)=(1-2*(bc==1)).*V1L(Imin+2);
        VL(Imax+3)=(1-2*(bc==1)).*V1L(Imax+2);
        %Update boundaries for Lax-Wendroff
        ULW(Imin+1)=ULW(Imin+2)+A*sin(2*pi*((n*dt)/(Xdomain/(C*M))))*exp(-.5*((n*dt-Tmax*T)^2)/(Tmax*G)^2);
        ULW(Imax+3)=ULW(Imax+2);
        VLW(Imin+1)=(1-2*(bc==1))*VLW(Imin+2);
        VLW(Imax+3)=(1-2*(bc==1))*VLW(Imax+2);
    end;
    
    %Show Results    
    plot(x,Ua,'k:',x,ULW,'k--',x,UL,'k-',x,Va,'b:',x,VLW,'b--',x,VL,'b-');
    axis([0 100 -4*A 4*A]);
    xlabel('Horizontal Position');
    ylabel('Amplitude');
    title(['Lax (Solid), Lax-Wendroff(Dashed) and Initial Condition (Dotted) @ t=',num2str(n*dt), ', n=',num2str(n),'']);
    pause(0.005); 
    
end;
