Solution (source code)

= Solution

For a normalized <wavefunction>, the <probability density> and <probability current> are
$$
\rho=|\psi|^2,\qquad \mathbf J=\frac{\hbar}{2mi}\left(\psi^*\nabla\psi-\psi\nabla\psi^*\right).
$$
The <Schrödinger equation> and its conjugate, for a real potential, are
$$
i\hbar\psi_t=-\frac{\hbar^2}{2m}\nabla^2\psi+V\psi,\qquad
-i\hbar\psi_t^*=-\frac{\hbar^2}{2m}\nabla^2\psi^*+V\psi^*.
$$
Multiplying respectively by $\psi^*$ and $\psi$ and adding their expressions for the time derivatives cancels the potential terms. Hence
$$
\rho_t=\frac{i\hbar}{2m}\left(\psi^*\nabla^2\psi-\psi\nabla^2\psi^*\right)=-\nabla\cdot\mathbf J.
$$
The gradient cross-terms cancel when expanding the divergence. Therefore
$$
\boxed{\rho_t+\nabla\cdot\mathbf J=0.}
$$
This <continuity equation> expresses local conservation of probability. Integrating over a region says that its probability decreases by the outward flux of the <probability current>.