L=300.;
N=2000;
Nt=300;
x=linspace(-L,L,N+1);
dx=2*L/N;
dt=5*dx;
a1=.6;
a2=.2;
v1=(1/a1-a1)/(1/a1+a1);
v2=(1/a2-a2)/(1/a2+a2);
g1=1/sqrt(1-v1^2);
g2=1/sqrt(1-v2^2);
x1=-L*.4;
x2=-L*.8;
phi12=zeros(Nt,N+1);
for k=1:Nt-2
phi12(k,:)=4*atan((a1+a2)/(a1-a2)*(exp(g1*(x-x1-v1*k*dt))...
    -exp(g2*(x-x2-v2*k*dt)))...
./(1+exp(g1*(x-x1-v1*k*dt)).*exp(g2*(x-x2-v2*k*dt))));
end
mesh(phi12)