clear all;

% Code to solve the KdV equation (with two coefficients a and b)
% u_t+a*uu_x+b*u_xxx=u_t+a/2*(u^2)_x+b*u_xxx=0
% based on code at:
% http://www.wikiwaves.org/Numerical_Solution_of_the_KdV
% with some modifications by A. Donev, Courant Institute

N = 256; % Number of grid points

if(1) % Start with two solitons and watch them collide
   % Here we are solving u_t+6*uu_x+u_xxx=u_t+3*(u^2)_x+u_xxx=0
   a=6; b=1;
   
   x = linspace(-10,10,N); % Domain
   c=16;
   u = 8*sech(2*(x+8)).^2 + 2*sech(x+1).^2;
   delta_t = 0.0001; % Time step size
   tmax = 1.0;
   limits = [-10 10 -1 10]; % axes limits for plot
   
else % Start with a Gaussian and watch solitons develop
  % See end of Maple worksheet NumericalSolutions.mw
  
  % Here we are solving u_t+uu_x+0.022^2*u_xxx=0
  a=1; b=0.022^2;
  
  x = linspace(-1,5,N);
  u = exp(-9*(x-1/2).^2);
  delta_t = 0.0005; % This is too large still for this method
  %delta_t = 0.0001; 
  tmax=6.0;
  limits = [-1 5 -0.25 2];

end

figure(1); clf
p = plot(x,u,'linewidth',3);
axis(limits);
grid on;

% ---------------------------------------------

delta_x = x(2) - x(1);
delta_k = 2*pi/(N*delta_x);

k = [0:delta_k:N/2*delta_k,-(N/2-1)*delta_k:delta_k:-delta_k];

nmax = round(tmax/delta_t);

U = fft(u);

for n = 1:nmax

   t = n*delta_t; % Current time
   
   % first we solve the linear part using Fourier series
   U = U.*exp(b*1i*k.^3*delta_t);

   %then we solve the non linear part using simple Euler method
   U = U - delta_t*((a/2)*1i*k.*fft(real(ifft(U)).^2));

   if mod(n,10)==0 % Plot every 10 steps
      u = real(ifft(U)); set(p,'ydata',u)
      title(sprintf('t = %7.5f',t),'fontsize',18), drawnow
   end   
   
end 
