I am trying to numerically solve the Gross-Pitaevskii equation for an impurity coupled with a one-dimensional weakly-interacting bosonic bath, given by (in dimensionless units):
\begin{align} i \frac{\partial \phi_0(x,t)}{\partial t} &= \left\{-\frac{1 + m}{2} \partial_{x}^2 + |\phi_0(x,t)|^2 + G\delta(x) + i V(t) \partial_{x} \right\}\phi_0(x,t)\nonumber\\ V(t) &= V_0 + i\frac{m}{\sqrt{\gamma}} \int dx \phi^*_0(x,t) \partial_{x} \phi_0(x,t) \end{align}
Where $m$,$\gamma$, $V_0$, and $G$ are external parameters which are real numbers. I would like to implement a conservative finite difference method, i.e., a finite difference method that discretizes the differential equation in both time $t$ and space $x$ and solve it while conserving the $L_2$ norm (number of particles), Energy and total momentum.
From studying different literatures (https://doi.org/10.3846/1392-6292.2009.14.109-126, https://www.mdpi.com/518624), I have found out that there is a finite discretization scheme (a ''modified version'' of the Crank-Nicholson scheme) that one can apply on a Gross-Pitaevskii equation which abides by the conservation rules. Without the drift term $i V(t)$, this scheme (derived using equation (3.1) of the first paper mentioned before) is given by:
\begin{align} i \frac{\phi_j^{n+1} - \phi_j^{n}}{\tau} = -\frac{1+m}{2} \left[ \frac{\phi_{j+1}^{n+1} - 2\phi_j^{n+1} + \phi_{j-1}^{n+1}}{2h^2} + \frac{\phi_{j+1}^n - 2\phi_j^n + \phi_{j-1}^n}{2h^2} \right] + \frac{1}{2}\left[\frac{G \delta(j)}{h} + \frac{1}{2} \left(|\phi_j^{n+1}|^2 + |\phi_j^n|^2 \right) \right] \left(\phi_j^{n+1} + \phi_j^n \right) \end{align}
Where $j$ and $n$ are discretized $x$ and $t$ respectively; and $h$ and $\tau$ denotes the values of space and time-steps, respectively. The initial conditions is given by $\phi^0_{j} = 1 \forall j \in [-L/2,L/2]$ and also periodic boundary conditions need to be implemented, given by $\phi_{L/2}^n = -\phi_{-L/2}^n \forall n \leq t_{\text{final}}$. I have two queries:
I have seen the Crank-Nicholson scheme to be implemented in the case of the heat-flow equation, where the linearity allows one to express the finite discretized equation to be written in the form of tridiagonal matrix and can be thus solved with Gaussian elimination processes (Thomas algorithm) (for a detailed example, check https://sepwww.stanford.edu/sep/prof/bei/fdm/paper_html/node15.html). But in the GPE case, the non-linearity prevents me from doing so. I would like to know how one can actually implement or solve this equation numerically.
Due to the physical nature of the problem, I would like the above-mentioned quantities to be also conserved in my case. With the addition of the drift term $iV(t)$, will the conservation laws still be respected with this conservative finite-difference scheme? I have checked that with the periodic boundary condition, the conservation laws are respected when there is no drift term.
Edit: I had a discussion with someone about these equations. Now I am curious if this set of non-linear coupled equations can be solved via Mathematica in some way?