N = 256
x = np.linspace(-10,10,N)
c=16
delta_x = x[2]- x[1]
delta_k = 2*np.pi/(N*delta_x);
k = np.arange(0,0.001, delta_k)
u = 1/2*c*(1/np.cosh(np.sqrt(c)/2*(x+8)))**2
delta_t = 0.4/N**2
tmax = 0.1
nmax = round(tmax/delta_t)
U = np.empty((len(x),nmax))
U[:,0] = np.fft.fft(u)
for i in range (1, nmax):
    #frst we solve the linear part
    U[:,i] = U[:,i-1]*np.exp(1j*k**3*delta_t);
    # then we solve the non linear part
    U[:,i] = U[:,i] - delta_t*(3j*k*np.fft.fft(np.real(np.fft.ifft(U[:,i]))**2))