spectral solution to 1D KdV equation

I am trying to solve the solution of a 1D KdV equation (ut+uux+uxxx=0) starting from a two solitons initial condition using Fourier spectral method. My code blows up the solution for a reason that I cannot figure out (probably due to a high order numerical stiffness!). I would appreciate any help. This is my main code file:

%Spatial variable on x direction
L=4; %domain on x
delta=0.05; %spatial step size
xmin=-L; %minimum boundary
xmax=L; %maximum boundary 
N=(xmax-xmin)/delta; %number of spatial points
x=linspace(xmin,xmax,N); %spatial vector
%--------------------------
%IC
sigma1 = 2; sigma2 = 4;
c1     = 2;     c2 = 1;
h1     = 1;     h2 = 0.3;
U = 0.*x;
U = h1*sech(sigma2*(x+c1)).^2+h2*sech(sigma1*(x-c2)).^2;
% %plotting
% figure(1)
% plot(x,U,'b','LineWidth',2);
% xlabel('$x$','Interpreter','latex')
% ylabel('$U(x,t=0)$','Interpreter','latex')
% xlim([-L-2 L+2])
% ylim([-0.5 1.2])
% set(gca,'TickLabelInterpreter','latex')
% title('Initial State')
% set(gca,'FontSize',16)
% drawnow;

%Fast Fourier Transform to the initial condition                 
Ut = fftshift(fft(U));   
Ut = reshape(Ut,N,1);
%--------------------------
%1D Wave vector disretisation              
k = (2*pi/L)*[0:(N/2-1) (-N/2):-1];
k(1) = 10^(-6);
k = fftshift(k);
k = reshape(k,N,1);

%first derivative (advection)
duhat = 1i .*k .*Ut;
du = real(ifft(ifftshift(duhat))); %inverse of FT

%third derivative (diffusion)
ddduhat = -1i .* (k.^3) .*Ut;
dddu = real(ifft(ifftshift(ddduhat))); %inverse of FT
%--------------------------
%Time variable
dt = 0.01; %time step
tmin = 0;
tmax = 4;
tspan = [tmin tmax];
%--------------------------
%solve
Time = 100; 
 for TimeIteration = 1:2:Time
    t= TimeIteration * dt;               
    %solve
    [Time,Sol] = ode113('FFT_rhs_1D',tspan,U,[],du,dddu); 
    Sol = Sol(TimeIteration,:);
    %plotting
    figure(2)
    plot(x,abs(Sol),'b','LineWidth',2);
    xlabel('$x$','Interpreter','latex')
    ylabel('$|{U(x,t)|}$','Interpreter','latex')
    xlim([-L-2 L+2])
    ylim([-0.5 1.2])
    set(gca,'TickLabelInterpreter','latex')
    set(gca,'FontSize',16)
    txt = {['t = ' num2str(t)]};
    text(1,1,txt,'FontSize',16)
    axis square
    drawnow
 end
 %--------------------------

And my function file is the follwoing:

 function rhs = FFT_rhs_1D(tspan,U,dummy,du,dddu)

%define diffusion/advection coeffieicnts
diffusion  = 1;
advection = 1;
%solve the right hand side
rhs = - advection .* U .* du - diffusion .* dddu;

end

Any help will be appreicated.

Thanks.

Lama

[Syed Abid Sahdman]

Did you get the solution?