% Spectral Method for Solving KdV 
% Modification of https://wikiwaves.org/Numerical_Solution_of_the_KdV
clear
N = 256; 
x = linspace(-10,10,N); 
delta_x = x(2) - x(1); 
delta_k = 2*pi/(N*delta_x); 
k = [0:delta_k:N/2*delta_k,-(N/2-1)*delta_k:delta_k:-delta_k]; 

% One soliton with speed c
c = 10; 
u = 1/2*c*(sech(sqrt(c)/2*(x+5))).^2; 

% Add other forms of initial conditions
% One soliton
eta0 = 1;
x0 = 5.0;
u = funckdv(x,0,eta0,x0);

% Two soliton
u = twosol(x,0);

delta_t = 0.4/N^2; 
tmax = 1; 
nmax = round(tmax/delta_t); 

nsteps=1000;
loops=floor(nmax/nsteps);
j=0;

U = fft(u);

for n = 1:nmax 
    % Solve the linear part 
    U = U.*exp(1i*k.^3*delta_t);
    % Solve the non linear part 
    U = U - delta_t*(3i*k.*fft(real(ifft(U)).^2)); 
    u1=real(ifft(U));
    if mod(n,nsteps)==0 
        j = j+1;
        plot(x,u1)
        axis([-10,10,-1,8])
        drawnow
        F(j) = getframe;        
    end 
end 

% Use below command to see the movie
% movie(F)
