function KdVrkm1()
global freq freq3
% solves u_t + u_{xxx} + (0.5u^2)_x = 0
tmax = 500; 
dt = 0.1; % step in time
n = 256;
x = -n : n; % mesh points
N = 2*n + 1; % the # of mesh points  = the length of the interval L
freq = -2i*pi*[-n : n]/N; % frequencies times i
u = zeros(1,N);
%% if you want to run  a test set flag = 1, else set flag = 2

flag = 2;

figure; clf; 
hold on; grid;

% initial data
if flag == 1
    u0 = 1./(cosh(x/sqrt(12))).^2; %                      (*)
    % plot of the exact solution for u_0(x) given by  (*)
    hp = plot(x,u0,'LineWidth',2); 
else
    u0 = exp(-(x/20).^2) +exp(-((x+100)/8).^2);
end
hpic = plot(x,u0,'LineWidth',2,'color','r'); % plot of the numerical solution
drawnow;

t = 0;
freq3 = freq.^3;
e3 = exp(freq3*dt);
em3 = exp(-freq3*dt);
while  t < tmax 
    t = t + dt;
    vk = fftshift(fft(u0)); % solve v_t = em3*N(e3*v) in the Fourier domain
    k1 = rhs(0,vk);
    k2 = rhs(0.5*dt,vk + 0.5*dt*k1);
    k3 = rhs(0.5*dt,vk + 0.5*dt*k2);
    k4 = rhs(dt,vk + dt*k3);
    vkn = vk + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
    un = ifft(ifftshift(e3.*vkn)); % return from the Fourier domain
    set(hpic,'xdata',x,'ydata',un);
    drawnow;
    u0 = un;
    if flag == 1
       set(hp,'xdata',x,'ydata',1./(cosh((x-t/3)/sqrt(12))).^2);
    end  
end
end
%%
function rk = rhs(dt,v);
global freq freq3
e3 = exp(freq3*dt);
em3 = exp(-freq3*dt);
vk1 = e3.*v;         %  exp(tL)v in Fourier space
v2 = ifft(ifftshift(vk1));      % exp(tL)v in regular space
v2 = 0.5*v2.^2;          % (exp(tL)v)^2/2 in regular space
rk = em3.*freq.*fftshift(fft(v2)); % exp(-tL)v(N(exp(tL)v)) in Fourier space
end




