%%
close all hidden
clear all
clc
set(gca,'FontSize',18)
set(gca,'LineWidth',2)

N = 8*256;
x = linspace(-40,40,N);
delta_x = x(2) - x(1);
delta_k = 2*pi/(N*delta_x);

k = [0:delta_k:(N/2-1)*delta_k,0,-(N/2-1)*delta_k:delta_k:-delta_k];

b = 20;epsilon = 1;
u = u0/(2*epsilon)*(- tanh(10*(x-epsilon)) + tanh(10*(x+epsilon)));

y_max = 40 ;

name = 'wide_function.gif';

eval(['delete ',name])

delta_t = 0.4/N^2;
t=0;
plot(x,u,'LineWidth',2)
axis([-10 10 -y_max/3 y_max])
xlabel('x')
ylabel('u')
text(5,0.7*y_max,['t = ',num2str(t,'%1.3f')],'FontSize',18)
drawnow
gif_add_frame(gcf,name,2);

tmax = 0.1; nplt = floor((tmax/20)/delta_t); nmax = round(tmax/delta_t);
udata = u.'; tdata = 0;
U = fft(u);

for n = 1:nmax
    t = n*delta_t;
    % first do the linear part
    U = U.*exp(1i*k.^3*delta_t);
    
    % then solve the nonlinear part
    
    U = U  - (3i*k*delta_t).*fft((real(ifft(U))).^2);
    
    if mod(n,nplt) == 0
        u = real(ifft(U));
        udata = [udata u.']; tdata = [tdata t];
        plot(x,u,'LineWidth',2)
        axis([-10 10 -y_max/3 y_max])
        xlabel('x')
        ylabel('u')
        text(5,0.7*y_max,['t = ',num2str(t,'%1.3f')],'FontSize',18)
        drawnow
        gif_add_frame(gcf,name,2);
    end
end

figure
%%

waterfall(x(769:1:1280),tdata(1:end),udata(769:1:1280,1:end)')
colormap(1e-6*[1 1 1]); view(-20,25)
xlabel x, ylabel t, 
% axis([-10 10 0 tmax y_max/3 y_max]), 
grid off
set(gca,'ztick',[0 y_max]), pbaspect([1 1 .13])

print -djpeg widefunction