clear

syms x t v1 v2 x1 x2
F = 4*atan(exp((x-x1-v1*t)/sqrt(1-v1^2))) ... 
    + 4*atan(exp(-(x-x2+v2*t)/sqrt(1-v2^2)));
%F = 4*atan(exp((x-x1-v1*t)/sqrt(1-v1^2)));
u = diff(F,x);
v = diff(F,t);

% Convert to MATLAB Functions @(t,v1,v2,x,x1,x2)

psi = matlabFunction(F);
U   = matlabFunction(u);
V   = matlabFunction(v);

% Parameters
x1 =-5;
x2 = 5;
v1 = 0.8;
v2 = 0.8;

M = 600;
xmin = -15;
xmax = 15;
h = (xmax - xmin)/M;
N = 300;
x = xmin:h:xmax;
% Note dt = dx = h;
T = N*h;

% Initial Conditions
psi0 = psi(0,v1,v2,x,x1,x2);
u0 = U(0,v1,v2,x,x1,x2);
v0 = V(0,v1,v2,x,x1,x2);
%psi0 = psi(0,v1,x,x1);
%u0 = U(0,v1,x,x1);
%v0 = V(0,v1,x,x1);


% Time Step
for i = 1:N
    [u1, v1, psi1]=sGstep(u0, v0, psi0, h);
    %[u2, v2, psi2]=sGstep(u1, v1, psi1, h);
    plot(x,psi1)
    axis([-15,15,0,15])
    pause(.25)
    u0 = u1;
    v0 = v1;
    psi0=psi1;
end
