% 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-(alpha*UXXX+beta*U.*UX)*(2*dt);

clear

% KdV Parameters
% Interval= [a,b]
% Tsteps  = Number of time steps 
% N       = Number of spatial subintervals
% delta   = dispersion coefficient
% beta    = nonlinear coefficient

a=0;
b=2;  
Tsteps=4585990; 
N=200;
N1=N+1;
delta = 0.022;
alpha = delta^2;
beta  = 1.0;

% Time and space increments
dx=(b-a)/N;
dt=dx^3/4;
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);

% Periodic initial condition
uex=cos(pi*x)'; 

% 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 - (alpha*(up2-2*up1+2*um1-um2)*dt/dx^3 ... 
        +  beta*(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 - (alpha*(up2-2*up1+2*um1-um2)*dt/dx^3 ... 
          +  beta*(up1+u0+um1)/3.*(up1-um1)*dt/dx);

  % Update u0 and u1 for next time step.
  u0=u1;
  u1=u2;
  if i== 1273240
      ubreak=u0;
  end
  if i== 4583661
      utrain=u0;
  end
  if mod(i,5000)==0 
    plot(x,u2)
    axis ([0,2,-1,3])
    title(['t = ', num2str(i*dt)])
    drawnow
  end
en(i)=sum(u2.^2);
end;

figure
plot(x,uex,'-.r', x,ubreak,'--k', x,utrain)
axis ([0,2,-1,3])
xlabel ('x')
ylabel('u(x,t)')
title(['Emergence of Solitons'])

% Embelish figure
hleg = legend('Initial wave','Shock formation','Emerged solitons',...
              'Location','NorthEast')
set(hleg,'FontAngle','italic','TextColor',[.3,.2,.1])
grid on
set(gca, 'GridLineStyle', '-');
grid(gca,'minor')