# Using Crank-Nicolson to solve Non-Linear Schrödinger equation in Python

I aim to solve the (non-linear) Schrodinger equation using the Crank-Nicolson method in Python. Here are my two functions.

def diags(N, dt, dx):
Asup = sp.empty(N)
Asub = sp.empty(N)
bdiag = sp.empty(N)
bsup = sp.empty(N)
bsub = sp.empty(N)

Asup.fill(-dt/(2*dx**2))
Asub.fill(-dt/(2*dx**2))
bdiag.fill(1 - dt/dx**2)
bsup.fill(dt/(2*dx**2))
bsub.fill(dt/(2*dx**2))

A = sp.sparse.spdiags([Adiag, Asup, Asub], [0, 1, -1],N, N)
b = sp.sparse.spdiags([bdiag, bsup, bsub], [0, 1, -1],N, N)
return A, b

def crankNicolson(psi, dt, V, A, b, g=0):
psi = np.exp(-dt*(V + g*np.abs(psi.reshape(-1))**2))*psi
psi = sp.sparse.linalg.bicg(A, b*psi.reshape(-1).astype(np.complex128))[0]
return psi


which has been taken from here. However, the given link aims to propagate in imaginary time, therefore I would need to have $$dt \rightarrow -i dt$$ in order to convert the code to propagate in real time instead. This can simply be done by putting -1j*dt as an input into the functions instead of just dt.

To put simply, the function diags() is run once to obtain A and b, which are then inputs for the function crankNicolson(), which in turn is run within a loop to iterate the process and propagate it in time.

However when I do this and propagate, my situation does not evolve in time but remains stationary. Can anyone see what's up with it? It absolutely should evolve, as I compare it to an explicit-midpoint method of mine which evolves in time well using the exact same input parameters.

The equation in dimensionless form is $$i\frac{\partial \psi}{\partial t} = -\frac{1}{2}\frac{\partial^2}{\partial x^2}\psi + V \psi + g|\psi|^2\psi$$
where the $$V$$ term is a self-consistent potential which is calculated separately (it is not important for this question as it is definitely correct), and $$g$$ is the coupling constant. The boundary conditions are that $$\psi(0) = 0$$ and $$\psi(L) = 0$$, owing to the fact that this is a decomposition of a spherically symmetric system with wavefunction $$\phi$$ according to $$\phi = \psi/x$$ (though this detail is perhaps not important).
The initial condition is a half-Gaussian, given by $$\psi = x \alpha \exp(-x^2/\beta)$$, which is centred on $$x=0$$.