clear all
close all
%------------------------------
% Solution of the KdV equation
%  u_t=6*u*u_x-u_xxx, -L<x<L,
%  Initial condition is
%  u(x,0)=-6*sech^2(x) and
% periodic boundary conditions, i.e., 
% u(-L,t)=u(L,t).L=1; % length of the string
%------------------------------
L=10
M=500;%600 % amount of points on the (-L,L)
x=linspace(-L, L, M+1);
dx=x(2)-x(1);% space grid size
x1=linspace(-L-2*dx, L+2*dx, M+5);
dx1=x1(2)-x1(1);
dt=2e-5; %5e-4;
T=100000;
tt=1000;
a=dt/dx;
%-----------------------------
u=zeros(T+1,M+5);
f=inline('-6*sech(y).^2', 'y'); % an initial condition for two-solitons solution u(x,0)=-N(N+1)sech^2(x), where N=2
u(1,:)=f(x1);
plot(x1,u(1,:),'r','LineWidth',2)
xlabel('x','FontName','Times New Roman','FontSize',16)
ylabel('u(x)','FontName','Times New Roman','FontSize',16)
axis([-L L -10 2])
pause(0.5)
%-----------------------------
u(1,1)=u(1,M+1);
u(1,2)=u(1,M+2);
u(1,3)=u(1,M+3);
u(1,M+4)=u(1,4);
u(1,M+5)=u(1,5);
%-----------------------------
 for i=3:M+3 % FTCS for the first time step
  u(2,i)=u(1,i)+3*a*u(1,i)*(u(1,i+1)-u(1,i-1))-a*(u(1,i+2)-2*u(1,i+1)+2*u(1,i-1)-u(1,i-2))/(2*dx^2);
 end;
 plot(x1,u(2,:),'r','LineWidth',2)
xlabel('x','FontName','Times New Roman','FontSize',16)
ylabel('u(x)','FontName','Times New Roman','FontSize',16)
axis([-L L -10 2])
u(2,1)=u(2,M+1);
u(2,2)=u(2,M+2);
u(2,3)=u(2,M+3);
u(2,M+4)=u(2,4);
u(2,M+5)=u(2,5);
%  %-----------------------------
for t=2:T % modifed time steping for t=2,3,....T
  for i=3:M+3
      u(t+1,i)=u(t-1,i)+2*a*(u(t,i+1)+u(t,i)+u(t,i-1))*(u(t,i+1)-u(t,i-1))-a*(u(t,i+2)-2*u(t,i+1)+2*u(t,i-1)-u(t,i-2))/(dx^2);
  end;
  if mod(t,tt)==0
       plot(x1,u(t+1,:),'r','LineWidth',2);
       axis([-L L -10 2])
       pause(0.5) 
  end;
u(t+1,1)=u(t+1,M+1);
u(t+1,2)=u(t+1,M+2);
u(t+1,3)=u(t+1,M+3);
u(t+1,M+4)=u(t+1,4);
u(t+1,M+5)=u(t+1,5);
end;
%--------------plot 3D plot--------------
tz=1:tt:T+1;
x2=x(3:M);
% figure(2)
% h=surf(x2,tz,u(tz,3:M));
%  shading interp
%  lighting phong
%  alpha(0.8)
figure(2)
[C,h]=contourf(x2,tz,u(tz,3:M),50);
set(h,'LineStyle','none')
colormap('jet')
xlabel('x','FontName','Times New Roman','FontSize',16)
ylabel('t','FontName','Times New Roman','FontSize',16)

