Solution (source code)

= Solution

Varying the potential gives
$$
\delta F_{\mu\nu}=\partial_\mu\delta A_\nu-
\partial_\nu\delta A_\mu.
$$
Antisymmetry implies
$$
\frac12F^{\mu\nu}\delta F_{\mu\nu}
=F^{\mu\nu}\partial_\mu\delta A_\nu.
$$
After integration by parts and omission of the boundary term,
$$
\delta S=-\frac1{\mu_0c}\int
\left[-\partial_\mu F^{\mu\nu}+m^2A^\nu-\mu_0J^\nu\right]
\delta A_\nu\,d^4x.
$$
Stationarity for arbitrary compactly supported $\delta A_\nu$ therefore gives the <Proca equation>
$$
\boxed{\partial_\mu F^{\mu\nu}-m^2A^\nu=-\mu_0J^\nu.}
$$