% Solves KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
% FFT with integrating factor phihat = exp(-ik^3t)*uhat.
% 
clear all
global k  k3

N = 256; dt = .4/N^2; x = (2*pi/N)*(-N/2:N/2-1)';
tmax=0.04;

A = 25; B = 16; clf, drawnow
%This is two soliton initial condition: a fat and thin one
u0 = 3*A^2*sech(.5*(A*(x+2))).^2 + 3*B^2*sech(.5*(B*(x+1))).^2;
%This is a Gaussian initial condition
%A=25;
%u0=A^2*exp(-x.^2/0.1);

%The state of the system at time (it-1)*dt is stored in
%as  ut(:,it) with it=1,2,....

%*****************************************
kx = [0:N/2-1 0 -N/2+1:-1]'; 
k=1i*kx;
k3=1i*kx.^3;
%k = 1i*kx; k3=1i*kx^3

uhat=fft(u0);uhat0=uhat;
phihat=uhat;
%in case we wanted to dealias
%nm=N/2+1;
%nmodes=N/2-1;nalias=floor(nmodes/3);
%dlia=ones(N,1);
%dlia(nm-nalias:nm+nalias)=zeros(2*nalias+1,1);
%dlia(nm)=1;
%phihat=phihat;


options = odeset('RelTol',1e-6,'AbsTol',1e-6); % Set tolerance
tspan=[0,dt];

Nt=floor(tmax/dt);
%Nt=5;
Tt=zeros(Nt,1);
ut=zeros(N,Nt);
ut(:,1)=u0;

for it=2:Nt;
    Tt(it)=(it-1)*dt;tspan=[(it-1)*dt,it*dt];
    
    [t,phihatt] =ode45('dtkdv', tspan, phihat, options); % Solve ODEs
    tn=length(t);
    phihat=phihatt(tn,:).';
    ut(:,it)=real(ifft(phihat.*exp(k3*Tt(it))));
    figure(3);
    plot(x,ut(:,it));%x,real(ifft(exp(-k*Tt(it)+k3*Tt(it)).*uhat0)),'o')
    title(['t = ',num2str(Tt(it))])
    drawnow
    %pause(0.1)
    
end;
gt=ut';
waterfall(x,Tt(1:50:Nt),gt(1:50:Nt,:)), colormap([0 0 0]), 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])




% set(gca,'ztick',[0 2000]), close(h), pbaspect([1 1 .13])