%
%  Solve 1D hyperbolic problem
%
%      u_t + a*u_x = b*u, 0<x<L, a>0
%      u(0,t)=0, t>0
%      u(x,0)=u0(x)
%
function lw(a,b,n,m)
if a <= 0
    fprintf(1,'Coefficient "a" should be positive\n');
    return;
end
L=2;
T=2;
dt=T/m;
x=L*(0:n-1)/(n-1);
h=L/(n-1);
c1 = dt*a/(2*h);
c2 = a*a*dt*dt/(2*h*h);
c3 = a*dt/h;
u=u0(x);
unew=zeros(n,1);
fprintf(1,'h=%.2f, dt=%.2f, Courant number = %.3f\n',h,dt,dt*(a-b*h)/h); 
for step=1:m
    %
    % Lax-Wendroff method
    %
    for j=2:n-1
        unew(j) = u(j) - c1*(u(j+1)-u(j-1)) + c2*(u(j+1)-2*u(j)+u(j-1)) + dt*b*u(j);
    end
    %
    % handle last node with backward differences:
    %
    unew(n) = u(n) - c3*(u(n)-u(n-1)) + dt*b*u(n);
    t = step*dt;
    uexact = exp(b*t)*u0(x-a*t); 
    if mod(step,1)==0
        plot(x,unew,x,uexact); axis([0 L -0.1 1]);
        str = sprintf('time = %f',t);
        title(str);
        pause(0.06);
    end
    u=unew;
end

%
% Initial value
%
function u=u0(x)
u=(x>=0.25 & x<=0.75).*cos(2*pi*(x-0.5)).^2;

