Solution (source code)

= Solution

Writing $D_A=\nabla+ieA/\hbar$, the current is
$$
J=-\frac{i\hbar}{2m}
(\psi^*D_A\psi-\psi(D_A\psi)^*).
$$
The two covariant factors acquire opposite phases under a gauge transformation, so $J$ is invariant.

Multiply the <Schrödinger equation> by $\psi^*$, subtract its complex conjugate multiplied by $\psi$, and use
$$
\psi^*D_A^2\psi-\psi(D_A^2\psi)^*
=\nabla\cdot\{\psi^*D_A\psi-\psi(D_A\psi)^*\}.
$$
The real scalar-potential terms cancel, leaving
$$
\partial_t|\psi|^2+\nabla\cdot J=0.
$$