% We generate solutions using the Zabusky-Kruskal scheme where
% the separate terms are given in terms of cyclically shifted vectors as
%   UXXX=(up2-2*up1+2*um1-um2)/(2*dx^3);
%   UX=(up1-um1)/(2*dx);
%   U=(up1+u+um1)/3;
% Then, ZK becomes
%   u2 = u0-(UXXX+6*U.*UX)*(2*dt);

clear

% KdV Parameters
% eta, x0 = soliton parameters
% Interval= [a,b]
% Tsteps  = Number of time steps 
% N       = Number of spatial subintervals

eta0=1;
x0=5.0;
a=-10;
b=10;
  
Tsteps=30000; 
N=200;
N1=N+1;

% Time and space increments
dx=(b-a)/N;
dt=dx^3/4/12;
t=(0:Tsteps-1)*dt;
x=linspace(a,b,N1);

% Array Initializations 
u0     = zeros(size(1:N1));
u1     = zeros(size(u0));
u2     = zeros(size(u0));
up2    = zeros(size(u0));
up1    = zeros(size(u0));
um1    = zeros(size(u0));
um2    = zeros(size(u0));
uex    = zeros(size(u0));
uexf   = zeros(size(u0));
en     = zeros(Tsteps,1);

% Select initial condition
    % One soliton solution
    uex=funckdv(x,0,eta0,x0)';
    uexf=funckdv(x,max(t),eta0,x0)';

    % Two soliton solution
    uex = twosol(x,0)'; 
    uexf = twosol(x,max(t))';

% Notation for different u's 
    % u0 = u at first time step.
    % u1 = u at second time step.
    % u2 = new value of u based on ZK scheme.

% First time step
% Compute cyclically shifted vectors
u0=uex;
um1=u0([ end 1:end-1 ]);
up1=u0([ 2:end 1 ]);
um2=u0([ end-1:end 1:end-2 ]);
up2=u0([ 3:end 1:2 ]);

% Apply ZK scheme
u1=u0-((up2-2*up1+2*um1-um2)*dt/dx^3+6*(up1+u0+um1)/3.*(up1-um1)*dt/dx)/2;

% Compute energy
en(1)=sum(u0.^2);
en(2)=sum(u1.^2);
              
   % Time Loop
for i=3:Tsteps 
  um1=u1([ end 1:end-1 ]);
  up1=u1([ 2:end 1 ]);
  um2=u1([ end-1:end 1:end-2 ]);
  up2=u1([ 3:end 1:2 ]);

  u2=u0-((up2-2*up1+2*um1-um2)*dt/dx^3+6*(up1+u0+um1)/3.*(up1-um1)*dt/dx);

  % Update u0 and u1 for next time step.
  u0=u1;
  u1=u2;
  if mod(i,100)==0 
    plot(x,u2)
    axis ([-10,10,-1,10])
    drawnow
  end
en(i)=sum(u2.^2);
end;

hold on
plot(x,uex,'r')
plot(x,uexf,'k')
hold off
