L=120.;
N=600;
Nt=1500;
x=linspace(-L,L,N+1);
dx=2*L/N;
dt=0.5*dx;
v=.35;
g=1/sqrt(1-v^2);
x0=L*.4;
phi=zeros(Nt,N+1);
phi(1,:)=tanh(g*(x-x0)/sqrt(2))-tanh(g*(x+x0)/sqrt(2))+1;
phi(2,:)=tanh(g*(x+v*dt-x0)/sqrt(2))-tanh(g*(x-v*dt+x0)/sqrt(2))+1;
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*(phi(k+1,:)-phi(k+1,:).^3);
end

Energy=zeros(Nt,N+1);
MaxIndex=zeros(Nt-1);
for i=1:Nt-1
Energy(i,:)=1/2*((phi(i+1,:)-phi(i,:))/dt).^2+...
    1/2*((phi(i,:)*shl(N)-phi(i,:))/dx).^2+1/4*(phi(i,:).^2-1).^2;
end
Energyhalf=zeros(Nt,N/2);
for j=1:N/2
    Energyhalf(:,j)=Energy(:,j);
end

for l=1:Nt-1
[C,MaxIndex(l)]=max(Energyhalf(l,:));
end

for m=1:Nt/3
Max1(m)=MaxIndex(m)*dx/dt;
Max2(m)=MaxIndex(Nt-m)*dx/dt;
end
%mesh(-phi)
mesh(Energy)
