Solution
= Solution
A backward <Euler method> with spatial <second-order central differences> gives
$$
\boxed{\frac{\Phi_j^{n+1}-\Phi_j^n}{\Delta\tau}=\frac{\Phi_{j+1}^{n+1}-2\Phi_j^{n+1}+\Phi_{j-1}^{n+1}}{h^2}-4\pi G\rho_j.}
$$
The source is fixed in pseudo-time; if it were evolving, its evaluation would also need specifying. Each pseudo-time step therefore requires a <linear system> solve rather than a purely local explicit update.