% 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 = .05; % 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.

  a =  2;
  b =  sqrt(a/2); % inverse of soliton width
  c =  -2*a; % soliton speed (for determining exact solution at final time)
  xm= 0.5*L; % Initial soliton center
  q = a*sech(b*(x-xm)).^2;

  tf = 1; % Final time
  nt = round(tf/dt); % Number of timesteps to take

  for it = 1:nt

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

    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

  Nf = 256;   % Number of plotting points
  xf = (0:(Nf-1))*L/Nf;
  qf = interpft(q,Nf);  % Numerical solution on plotting grid
  q0f = a*sech(b*(xf-xm)).^2;  % Exact initial soln on plotting grid
  qef = a*sech(b*(xf-xm-c*tf)).^2;  % Exact t=1 soln on plotting grid
  plot([x L],[q q(1)],'r+',[xf L],[qf qf(1)],'r-',...
       [xf L],[qef qef(1)],'b-',[xf L],[q0f q0f(1)],'b--') 
  xlabel('x')
  ylabel('q')
  title(['KdV soliton at  t = ' num2str(tf) 'using PS/RK4'])
  legend('Gridpoints','PS','Exact','Initial')