%one or two solitons colliding in various potentials. x-t plot
%variations: 
%sec = 0, k = 0, potent = 1 => off-centered solitons in parabolic potential
%sec = 1, k = 0, potent = 1 => watch two solitons attract in neg. potential
%sec = 1, k = 10, potent = 0,1,2,3 => watch soliton collisions in various potentials

% parameters ..........
n = 400;       %lattice points
dt = .00001;   %time step
potent = 0;    %potential type: 0 none, 1 parabolic, 2 step, 3 multi parabolic
T = .3;        %time duration
sec = 1;       %second wave. 0 or 1
k = 10;        %wave number

xmin = -6;
xmax = 6;
c = 1;          %width
p = 4;          %number of step parabolic potentials
x0 = 3;        %off centering of waves
s = 10;      %potential strength.
B = -2;        %non-linear strength. 
% .................

dx = (xmax-xmin)/n;
x = xmin + (1:n)'*dx;  % x is col vec

V = zeros(1,n);

if potent == 0;                 %no potential. potent = 0
else
   if potent == 1;
      V = s*x.^2;               %parabolic potential. potent = 1
   else
      if potent == 2;           %step potential. potent = 2
         for a = 3*n/8:5*n/8
            x0 = 4;
            V(a) = s;
         end
      else                      %multi parabolic potential. potent = 3
      V = s*multiParabolas(n, p, xmin, xmax, dx);
      end
   end
end

a = 1/(dx)^2 - V/2;
b = 1/2*1/(dx)^2;

% creation of K the sparse matrix
D = sparse(1:n,1:n,a,n,n);
E = sparse(2:n,1:n-1,b*ones(1,n-1),n,n);
E(1,n) = b;
K = E+D+E';

U = 2*sech(2*c*(x-x0)).*exp(-1i*k*x) + sec*2*sech(-2*c*(x+x0)).*exp(1i*k*x);
I = speye(n,n);           % keeps everything sparse 
t = 0;

%using CN scheme in paper...
Uold = U;
Unew = zeros(n,1);
modT = 100;
xtU = zeros(n,(T + modT*dt)/(modT*dt));     %records U as time goes on

while t < T
  if (mod(floor(t/dt),modT)==0)
    xtU(:,floor(t/(modT*dt) + 1)) = U.*conj(U);
    plot(x, [U.*conj(U)], x , V);  %real(U) imag(U)], x, V);
    axis([xmin xmax -1 10]);
    text(1,1, sprintf('t=%g, prob=%g', t, dx*dot(U,U)));
    drawnow;
  end
    Unew = (I - 1i*dt*K)\((I + 1i*dt*K)*U - B*1i*dt*(3/2*conj(U).*U.*U - 1/2*conj(Uold).*Uold.*Uold));
    Uold = U;
    U = Unew;
    t = t + dt;
end

figure;
surf(xtU);