%%
% Solve u_t+a*u_x=0 on [-inf,inf] with 3 schemes and a Rectangle IC
%%
clc;
clear all;
close all;
scheme='Lax-Wendroff'; %Upwind,Lax-Friedrichs,Lax-Wendroff
H=@ heaviside;
a = 1;T=1; ax=-5;bx=5;
errvals=[]; hvals=[];
for j=[100 200 400]
    
    h = (bx-ax)/j;
    k = 0.5*h;
    n = j+1;
    x=linspace(ax,bx,j+1);
      
    uo= 1 + H(x+1) - H(x-1);%initial condition with compact support
    u=uo;
    for t = 0:k:T
        I=2:n-1;
        switch(scheme)
            case 'Upwind'
                %Upwind
                u(I)=uo(I) - a*k*(uo(I)-uo(I-1))/h; %only update interior nodes
            case 'Lax-Friedrichs'
                %Lax-Friedrichs:
                u(I) = 0.5*(uo(I+1)+uo(I-1)) - a*k*(uo(I+1)-uo(I-1))/(2*h);
            case 'Lax-Wendroff'
                %Lax-Wendroff:
                u(I)=uo(I) - a*k*(uo(I+1)-uo(I-1))/(2*h)+...
                    (a*k)^2/(2*h^2)*(uo(I+1)-2*uo(I)+uo(I-1));
            otherwise
                error('Not a valid scheme!');
        end
        uo = u;
    %exact solution u(x,t)=eta(x-a*t)
    uexact = 1 + H(x+1-a*t) - H(x-1-a*t); %exact solution at time t 
    %plot each time step    
    plot(x,u,'--rx','LineWidth',2)
    hold on;
    plot(x,uexact,'-b','LineWidth',2)
    axis([ax bx 0 2.5])
    set(gca,'FontSize',20);
    xlabel('x','FontSize',20,'fontweight','b');
    ylabel('u','FontSize',20,'fontweight','b');
    legend(scheme,'Exact','Location','south');
    title(sprintf('mesh size [%d]',j));
    hold off
    
    pause(0.1);
    end
    %compute error at final time
    err = h*sum(abs(u-uexact)); %mean error, for discontinous solution
    
    errvals=[errvals;err];
    hvals=[hvals;h];
    
    pause(1) %next mesh size
    
end
fprintf('--------The order of accuracy: [%s] scheme-----------\n',scheme);
error_table(hvals, errvals);   % print tables of errors and ratios
