% Numerically calculate soln. to KdV eqn. q_t + 6qq_x +q_xxx= 0 on
% domain 0<x<1 with periodic BCs using spectral method with RK4
% time-differencing scheme. 
% Note that with RK4 the timestep-limiting term is
% the third derivative; since the maximum wavenumber is kmax = pi*(N/L), the 
% corresponding linear oscillation frequency for this waveno. is om = kmax^3 
% and the max stable frquency for RK4 is om*dt < 2.82, we have stability 
% restriction dt < cmax*(L/N)^3, where cmax = 2.82/pi^3 = 0.091.

  N = 32; % Number of Fourier modes
  L = 16;  % Domain size
  C = .002; % Nondimensional timestep parameter(gives dt = 0.05 for L/N=1)

  dt = C*(L/N)^3;
  x = L*(0:(N-1))/N;  %  x-gridpoints [1xN]
  dx = L/N;
  M = [0:(N/2-1) (-N/2):(-1)];  
  k = 2*pi*M/L;  % Wavenumbers [1xN].  

  % Specify soliton ICs. Amplitude a is user-defined. Note that a >> 1 to
  % ensure soliton width <<1, so periodic BCs don't substantially perturb the
  % soliton soln derived for an unbounded domain.

  a1 =  2;
  b1 =  sqrt(a1/2); % inverse of soliton width
%  c1 =  -2*a1; % soliton speed (for determining exact solution at final time)
  x1= 8; % Initial soliton center
  a2 =  4;
  b2 =  sqrt(a2/2); % inverse of soliton width
%  c2 =  -2*a2; % soliton speed (for determining exact solution at final time)
  x2= 12; % Initial soliton center
  q = a1*sech(b1*(x-x1)).^2 + a2*sech(b2*(x-x2)).^2 % Soliton IC.

  % Plot initial condition for x-t offset plot

  yscale = 0.02;
  tp = 0.1; % Time between plots
  tf = 1.5; % Final time

  Nf = 256;   % Number of plotting points
  xf = (0:(Nf-1))*L/Nf;
  qf = interpft(q,Nf);  % Numerical solution on plotting grid
  plot(xf,yscale*qf);
  xlabel('x')
  ylabel(['t + ' num2str(yscale) 'q(x,t)'])
  title(['Interacting KdV solitons, PS, N=' num2str(N)])
  text(0.05*L,1.05*tf,['dt=' num2str(dt)])
  axis([0 L 0 1.1*tf])
  hold on

  np = round(tf/tp); % Number of times to plot
  nt = round(tp/dt); % Number of timesteps between plots

  for ip = 1:np
    for it = 1:nt

      % March forward dq/dt = -S_KdV(q) using RK4
      % S_KdV(q) = -6qq_x - q_xxx is evaluated pseudospectrally.

      d1 = -dt*S_KdV(q,L); 
      d2 = -dt*S_KdV(q + 0.5*d1,L);
      d3 = -dt*S_KdV(q + 0.5*d2,L);
      d4 = -dt*S_KdV(q + d3,L);
      q = q + (d1 + 2*d2 + 2*d3 + d4)/6; % q marched forward dt
    end
    qf = interpft(q,Nf);  % Numerical solution on plotting grid
    plot(xf,yscale*qf+ip*tp);
  end
  hold off