% p1c.m - Solve KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(-ik^3t)*u-hat.
%         u(x,0) = 3A^2sech^2(A(x+2)/2)+3B^2sech^2(B(x+1)/2),A=25,B=16

% Set up grid and two-soliton initial data:
  N = 256;                           % number of spatial grid points
  h = 2*pi/N;                        % spatial grid size
  x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
  dt = 25/N^3;                       % stable time step
  dt = 30/N^3;                       % unstable time step
  A = 25; B = 16;                    % parameters for initial value
  u = 3*A^2*sech(.5*(A*(x+2))).^2 + ...
      3*B^2*sech(.5*(B*(x+1))).^2;   % initial value of u
  v = fft(u);
  k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
  ikv = (1i*k).^3;                   % scaling factor for 3rd partial deriv.

% Solve PDE and plot results:
  tmax = 0.006;                      % end time
  nplt = floor((tmax/25)/dt);        % number of solutions to plot
  nmax = round(tmax/dt);             % number of solutions to compute
  udata = u;                         % keep track of solutions
  tdata = 0;                         % keep track of time steps
  clf, drawnow;                      % clear plot
  set(gcf,'renderer','zbuffer');     % set rendering more for fast plotting
  H = gcf;                           % get handle of current plot
  h = waitbar(0,'please wait...');   % create waitbar
  t = 0;
  E = exp(t*ikv);                    % integrating factor at t
  Edt2 = exp(-dt*ikv/2);             % integrating factor at dt/2 
  Edt = exp(-dt*ikv);                % integrating factor at dt
  for n = 1:nmax                     
     t = n*dt;                       % set time for solution 
     g = -1;                         % scaling factor of RHS 

     % 4th-order Runge-Kutta method
     F1 = dt*g.*fft(f1c(v,k));
     F2 = dt*g.*Edt2.*fft(f1c(v+F1/2,k));
     F3 = dt*g.*Edt2.*fft(f1c(v+F2/2,k));
     F4 = dt*g.*Edt.*fft(f1c(v+F3,k));
     v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;

     if mod(n,nplt) == 0             % plot the solution 
      u = real(ifft(v));             % compute solution in real domain
      waitbar(n/nmax);               % update waitbar
      udata = [udata u];             % keep track of plotted data
      tdata = [tdata t];             % keep track of plotted time
    end
  end
  
  % Trefethen's code for plotting
  waterfall(x,tdata,udata'), colormap(1e-6*[1 1 1]); view(-20,25)
  xlabel x, ylabel t, axis([-pi pi 0 tmax 0 2000]), grid off
  set(gca,'ztick',[0 2000]), close(h), pbaspect([1 1 .13]),figure(H);