clear all
close all
%-----------------------------------------------------------------------
% Solve KdV-Equation u_t-6uu_x+u_xxx=0 on [-pi, pi] by pseodo-spectral
% method with integrating factor exp(-i*k^3*t)
%
%------------------------------------------------------------------------

% set the grid and a two-soliton initial data
N=256; 
dt=0.4/N^2; 
x=(2*pi/N)*(-N/2:N/2-1)';
A=25; B=16;
u=A^2*sech(0.5*(A*(x+2))).^2/2+B^2*sech(0.5*(B*(x+1))).^2/2; %initialcondition
plot(x,u,'r','LineWidth',2)
xlabel x, ylabel u
axis([-pi pi 0 310])
%-----------------
v=fft(u); 
k=[0:N/2-1 0 -N/2+1:-1]';
ik3=1i*k.^3;
%--------------------------------------------------------------------
tmax=0.006;
nplt=floor((tmax/25)/dt);
nmax=round(tmax/dt);
%----------------
udata=u; tdata=0; 
h=waitbar(0,'please wait...');
%--------------------------------------------------------------------
% Time-Loop with RK4
for n=1:nmax
    t=n*dt;
    g=-3i*dt*k;
    E=exp(dt*ik3/2);
    E2=E.^2;
    a=g.*fft(real(ifft(v)).^2);
    b=g.*fft(real(ifft(E.*(v+a/2))).^2);
    c=g.*fft(real(ifft(E.*v+b/2)).^2);
    d=g.*fft(real(ifft(E2.*v+E.*c)).^2);
    v=E2.*v+(E2.*a+2*E.*(b+c)+d)/6;
    if mod(n,nplt)==0
        u=real(ifft(v));
        waitbar(n/nmax);
        udata=[udata u];
        tdata=[tdata t];
    end;
end;
%---------------------------------------------------------------------
% Plot results
waterfall(x,tdata,udata'),
colormap([1 0 0]),
view(-20,25)
xlabel x, ylabel t, axis([-pi+0.01 pi-0.01 0 tmax 0 500]), grid off
close(h), pbaspect([1 1 0.13])
