Currently, I am studying a 2D nonlinear Schroedinger equation and searching for the fastest method. $$ \begin{equation} i \frac{\partial \psi}{\partial t} = [-\frac{1}{2} \nabla^2 + V_0(r) - i\gamma + g|\psi|^2]\psi + \Phi(r,t), \end{equation} $$ Studying the amplitude of the solution to this time-dependent equation I decided to split the $\psi$ function into two independent functions of time and coordinate: $$ \begin{equation} \psi = u(r)e^{-i\varepsilon t}, \Phi(r,t) = A(r)e^{-i\varepsilon t} \end{equation} $$ This simplification speeds up the solution program because I no longer need to simulate the evolution of the system and just observe the final result at the stationary regime. Then I obtain such equation: $$ \begin{equation} [-\frac{1}{2} \nabla^2 + V_0(r) - i\gamma - \varepsilon + g|\psi|^2]\psi + A(r) = 0 \end{equation} $$
I use Newton's method to find the solution (note that $r$ stands for $(x,y)$). Because of the 2D problem, I discretize the domain in 1024x1024 vertices. This leads to a problem of inversing (10241024)x(10241024) matrix. This task doesn't allow me to use GPU because I have low memory resources and I can't parallelise it with CPU.
At this moment a MATLAB code I have created takes way too much time (up to $2500$ seconds to converge to accuracy $10^{-9}$). So I wonder if there is a faster method I could implement to achieve the same accuracy within a smaller time.
