Hi I am trying to solve a variation of the heat equation with interaction terms and external source — actually it’s the Schrodinger-Newton equation if that’s more familiar to you. \begin{align} i \dot \psi & = -\nabla^2 \psi + \psi \phi_N + |\psi|^2\psi, \cr \nabla^2 \phi_N & = |\psi|^2, \end{align} The boundary condition is that the $ \psi(r)$ vanishes at large radius and its first derivative vanishes at zero. The initial condition is not known and the only requirement is the wave function being real at t=0.

Since NDSolve requires an explicit initial condition, I’ve been trying the following code with a guessed initial condition $ \psi \sim \mathrm{e}^{-r}$ but without much success. The kernel complains about inconsistency (ibcinc), insufficient number of boundary conditions (bcart), unable to continue with complex values (mconly), and unable to find initial conditions that satisfy (icfail). Is there any way to solve it with NDSolve?

tmp\[Psi]initial[r_] := (1 + r)*E^(-r) tmp\[Phi]initial[rr_] := (2/rr*D[tmp\[Psi]initial[r], r] + D[tmp\[Psi]initial[r], r, r]) /. r -> rr rStart0 = 0.01; rEnd0 = 5;  eqn = {-I*D[\[Psi][r, t], t] - 1/2*(D[\[Psi][r, t], r, r]          + 2/r*D[\[Psi][r, t], r]) + \[Psi][r, t]*\[Psi][r, t]* Conjugate[\[Psi][r, t]] + \[Psi][r, t]*\[Phi][r, t] == 0 , 2/r*D[\[Phi][r, t], r] + D[\[Phi][r, t], r, r] == Conjugate[\[Psi][r, t]]*\[Psi][r, t]};  ic = {\[Psi][r, 0] == tmp\[Psi]initial[r], \[Phi][r, 0] == tmp\[Phi]initial[r]}; bc = {(Derivative[1, 0][\[Psi]][rStart0, t]) == (D[tmp\[Psi]initial[r], r] /. r -> rStart0)  , (\[Psi][r, t] /. r -> rEnd0) == tmp\[Psi]initial[rStart0]  , \[Phi][rStart0, t] == tmp\[Phi]initial[rStart0] , \[Phi][rEnd0, t] == tmp\[Phi]initial[rEnd0]};  sol = NDSolveValue[Flatten@{eqn, ic, bc}, {\[Psi], \[Phi]}, {r, rStart0, rEnd0}, {t, 0, 10}, MaxSteps -> 10^6]