Solution (source code)

= Solution

For the supplied <Lagrangian>, variation with respect to $\phi^*$ uses
$$
\frac{\partial\mathcal L}{\partial(\partial_\mu\phi^*)}
=\partial^\mu\phi+ieA^\mu\phi,
\qquad
\frac{\partial\mathcal L}{\partial\phi^*}
=-m^2\phi-ieA_\mu\partial^\mu\phi+f.
$$
Thus $\Box\phi+2ieA^\mu\partial_\mu\phi+ie(\partial_\mu A^\mu)\phi+m^2\phi=f$. In the <Lorenz gauge>,
$$
\boxed{(\Box+2ieA^\mu\partial_\mu+m^2)\phi=f(t).}
$$
For real external fields and driving, the conjugate equation has $-2ieA^\mu\partial_\mu$ acting on $\phi^*$. The supplied model has no $e^2A^2|\phi|^2$ term; one must not insert such a term from a different minimally coupled <Lagrangian>.