nonlinear Schroedinger equation resolution with split step method

Question

I'm trying to resolve a nonlinear PDE Schrödinger equation with the split-step Fourier pseudo-spectral method. The fonction "solver" is supposed to solve the time and space dependant equation with the split-step method but it gives values oscillating between 0.49 and 0.51 instead of the good answer which is between 0 and 1.657 (see here). What should I change to get the good result?

import numpy as np
from matplotlib import pyplot as plt


i = np.complex(0,1)

K = 1

dt = 0.01 #temps
t = np.arange(0,20,dt)

N = 1024#espace
x = np.linspace(-20,20,N) #domaine truncated
L = 40
dx = L/N

k =np.linspace(-N/2,N/2-1,N)

def solver(x,t):
    F = []
    psi = 0.5 +0.01*np.cos(2*np.pi*x/40)
    for l in range(len(t)):
        g = np.exp(-i*K*(np.abs(psi))*dt)*psi
        fg = np.fft.fftshift((np.fft.fft(g)))
        fpsi = np.exp(-i/2*(2*np.pi*k/L)**2*dt)*fg
        psi = np.fft.ifft(fpsi)
        F.append(psi)
    return F

M = np.abs(solver(x,t))
print(M)
#plt.plot(x,M)
#plt.show()

fig = plt.figure(figsize=(10,10))
ax1 = fig.add_subplot()
m1 = ax1.matshow(M)
plt.title("")
fig.colorbar(m1)
ax1.set_xlabel('espace')
ax1.set_ylabel('temps')

enter image description here