L=120.;
N=600;
Nt=2000;
x=linspace(-L,L,N+1);
dx=2*L/N;
dt=0.5*dx;
v=.4;
g=1/sqrt(1-v^2);
x0=L*.4;
phi0=pi/5;
x1=-log(tan(phi0/4));
phi=zeros(Nt,N+1);
phi(1,:)=4*atan(exp(g*(-x-x0)))+4*atan(exp(x-x1));
phi(2,:)=4*atan(exp(g*(-x+v*dt-x0)))+4*atan(exp(x-x1));
for k=1:Nt-2
    phi(k+2,:)=2*phi(k+1,:)-phi(k,:)+dt^2/dx^2*(phi(k+1,:)*shl(N)...
        -2*phi(k+1,:)+phi(k+1,:)*shr(N))-dt^2*sin(phi(k+1,:));
   phi(k+2,N/2)=phi0;
end
mesh(phi)