% one-way wave equation u_t+au_x=0 a>0
%
a=2;
c=1;
% initial function
x=-30:0.25:30;
U = c./(2*cosh(0.5*sqrt(c).*x).^2)';
  plot(x,U)
  
N=100;
h=0.25;
k=0.124;
p=k/h;
%
mm=size(U,1);
S=sparse(diag(ones(mm-1,1),-1));
for n=1:N
    U=(1-a*p)*U+a*p*S*U;
    Uex = c./(2*cosh(0.5*sqrt(c).*(x-a*k*n)).^2);
    plot(x,U,'k')
    hold on
    plot(x,Uex,'r .')
    drawnow
    hold off
end

pause
U = c./(2*cosh(0.5*sqrt(c).*x).^2)';

for n=1:N
   
    Uex = c./(2*cosh(0.5*sqrt(c).*(x-a*k*n)).^2);
   
    plot(x,Uex,'r .')
    drawnow
    hold off
end

pause

U = c./(2*cosh(0.5*sqrt(c).*x).^2)';


Sc=sparse(-diag(ones(mm-1,1),-1)+diag(ones(mm-1,1),+1));
for n=1:N
    U=U-a*p/2*Sc*U;
    Uex = c./(2*cosh(0.5*sqrt(c).*(x-a*k*n)).^2);
    plot(x,U)
    hold on
    plot(x,Uex,'r .')
    drawnow
    hold off
end

pause
U = c./(2*cosh(0.5*sqrt(c).*x).^2)';

Sdd=diag(ones(mm-1,1),-1)+diag(ones(mm-1,1),+1);
A=sparse(-2*eye(mm)+Sdd);
for n=1:N
    U=U-a*p/2*Sc*U+1/2*(a*p)^2*A*U;
    Uex = c./(2*cosh(0.5*sqrt(c).*(x-a*k*n)).^2);
    plot(x,U,'g')
    hold on
    plot(x,Uex,'r .')
    drawnow
    hold off
end