# Boundary conditions for a Non-linear Schrödinger equation using an extended crank nicolson scheme

I try to solve numerically the following PDE for $$E(r, z)$$ with a cylindrical symmetrie (i. e. $$E(r, z) = E(-r, z)$$).

$$\frac{\partial E}{\partial z} = \frac{i}{2k} \Delta E + \mathcal{N}(E)$$

Where $$\Delta$$ is the Laplace operator in transversal direction and $$k$$ a real number. I want to use the Crank-Nicolson scheme to treat the Laplace operator and the Adams-Bashforth scheme to treat the nonlinearity ($$\mathcal N(E)$$). $$r$$ is defined by $$r_j = 0 + j\Delta r$$ for $$j = 0 \dots N$$. Since the Laplace in cylinder coordinates is given by: $$\Delta E = \frac{\partial^2E}{\partial r^2} + \frac{1}{r}\frac{\partial E}{\partial r}$$ one gets the following representation of the Laplace operator:

$$\Delta E_j^n = E^n_{j-1} -2 E^n_j + E^n_{j+1} + \frac{1}{2j}(E^n_{j+1} - E^n_{j-1})$$

The given PDE equation is therefore given by the following, where the nonlinearity is treated by the Adams-Bashforth scheme.

$$E_j^{n+1} - E_j^n = i\delta (\Delta_j E^{n+1}_j + \Delta_j E^{n}_j) + (3/2 \mathcal N^n_j - 1/2 N^{n-1}_{j})$$

where $$\delta = \frac{\Delta z}{4 k \Delta r^2}$$. From this $$E^{n+1}_j$$ can be expressed the following way:

$$E^{n+1}_j = L_-^{-1} [L_+E^n_j + 3/2\mathcal N^n_j - 1/2\mathcal N^{n-1}_j]$$ with the following matrixes $$L_-$$ and $$L_+$$.

$$L_{\pm} = \left( \begin{array}{rrrr} 1\mp2i\delta & \pm i\delta v_0 \\ \pm i\delta u_1 & 1\mp2i\delta & \pm i\delta v_1\\ & & & & \\ & & & & \\ & & & \pm i\delta u_N & 1\mp2i\delta\\ \end{array}\right)$$

Where $$u_j = 1 - 1/(2j)$$ and $$v_j = 1 + 1/(2j)$$.

For my problem the following boundary conditions are given: $$E(r = r_{max}, z) = 0$$ and $$\frac{\partial E}{\partial r} |_{r=0}$$. The first one for $$r=r_{max}$$ I can easily incorporate by changing the last row of $$L_{\pm}$$.

The other one ($$\frac{\partial E}{\partial r} |_{r=0}$$) gives me trouble. I somehow have to change the entries of the first row of $$L_{\pm}$$, but I don't know how. I got some working boundary conditions for the case the nonlinearity is 0, but they break as soon as I add a nonlinearity.

I am greatful for any help.