% kdv.m - solution of the Korteweg-DeVries equation by ETDRK4 scheme
%
% u_t = - u_xxx - uu_x with periodic boundary conditions on [-pi,pi]
% computation is based on v = fft(u), so linear term is diagonal
% see p27.m in Trefethen "Spectral Methods in Matlab"
%
% These are the functions as calculated after taking the FFT
% Linear Function L = i*k^3 (see notation in Kassam & Trefethen)
% Nonlinear Function N = -0.5*i*k(F(Re(F^-1(uhat))^2))

% Spatial grid and initial conditions
N = 256; 
x = (2*pi/N)*(-N/2:N/2-1)';
u = 3*25^2*sech(25.*(x+2)./2).^2 + 3*16^2*sech(16.*(x+1)./2).^2;
v = fft(u);

%Precompute various ETDRK4 scalar quantities
h = .4/N^2;    %fixed time step
k = [0:N/2-1 0 -N/2+1:-1]';      %wave numbers
L = 1i*(k.^3);
E = exp(h*L); E2 = exp(h*L/2);
M = 16;                             % no. of points for complex means
r = exp(1i*pi*((1:M)-.5)/M);        % roots of unity
LR = h*L(:,ones(M,1)) + r(ones(N,1), :);
Q = h*real(mean(    (exp(LR/2)-1)./LR    , 2));
f1 = h*real(mean( (-4-LR+ exp(LR).*(4-3*LR+LR.^2))./LR.^3   ,2));
f2 = h*real(mean( (2+LR+ exp(LR).*(-2+LR))./LR.^3   ,2));
f3 = h*real(mean( (-4-3*LR-LR.^2+ exp(LR).*(4-LR))./LR.^3   ,2));

% Main Time-stepping loop
uu = u; tt=0;
tmax = .006; nmax = round(tmax/h); nplt = floor((tmax/100)/h);

g = -0.5*i*k;
for n = 1:nmax
    t = n*h;
    Nv = g.*fft(real(ifft(v)).^2);
    a = E2.*v + Q.*Nv;
    Na = g.*fft(real(ifft(a)).^2);
    b = E2.*v + Q.*Na;
    Nb = g.*fft(real(ifft(b)).^2);
    c = E2.*a + Q.*(2*Nb-Nv);
    Nc = g.*fft(real(ifft(c)).^2);
    v = E.*v + Nv.*f1 + 2*(Na+Nb).*f2 + Nc.*f3;
    if mod(n,nplt)==0
        u = real(ifft(v));
        uu = [uu,u]; tt = [tt,t];
    end
end
%plotting code is borrowed from trefethen p27.m
waterfall(x,tt,uu'), 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]), pbaspect([1 1 .13])