Wat you want: Obtaining all four Maxwell equations from a Lagrangian in $F$, works perfectly well if you use the bi-spinor generators to define $F$ and the other electromagnetic fields instead. Using the chiral gamma matrices and the bi-spinor boost and rotation generators we get:
\[\begin{equation} \begin{array}{lrcl} \mbox{mass dimension 1:}~~~~ & \mathbf{A} &=& \gamma^\mu A^\mu \\ \mbox{mass dimension 2:}~~~~ & \mathbf{F} &=& \vec{K}\cdot\vec{E}-\vec{J}\cdot\vec{B} \\ \mbox{mass dimension 3:}~~~~ & \mathbf{J}\, &=& \gamma^\mu\,j^\mu \\ \end{array} \end{equation}\]
We can now write the laws of the Electromagnetic field, going from one mass dimension to another, by applying the differential operator matrix $/\!\!\!{\partial} =\gamma^\mu\partial_\mu = \sqrt{\Box}$ on the operator field matrices defined above.
\[\begin{equation} /\!\!\!{\partial}\mathbf{A}^{\!\dagger} = \mathbf{F}~~~~~~ ~~~/\!\!\!{\partial}\mathbf{F} = \mathbf{J}^\dagger \end{equation}\]
The complex conjugate transpose $\mathbf{A}^{\!\dagger}=\gamma^\mu A_\mu $ is the covariant form of $\mathbf{A}$.
In the first step we have applied the conservation law $\partial_\mu A^\mu\!=\!0$ on the diagonal and in the second step we find all four of Maxwell's laws, the inhomogeneous $\partial_\mu F^{\mu\nu}\!=\!j^\nu$ as well as the homogeneous $~\partial_\mu\! *\!\!F^{\mu\nu}\!=j^\nu_{^A}=\!0$.
We can write for the Lagrangian of the Electromagnetic field in vacuum:
\[\begin{equation} \mathcal{L} ~=~ \tfrac12\mathbf{F}\mathbf{F} \end{equation}\]
Where $\mathcal{L}$ is a matrix operator field invariant under Lorentz transform. To find the equations of motions: The four Maxwell equations. We write:
\[\begin{equation} \mathcal{L} ~=~ \tfrac12/\!\!\!{\partial}\mathbf{A}^{\!\dagger}/\!\!\!{\partial}\mathbf{A}^{\!\dagger}~~~~ \underrightarrow{\mbox{ Euler Lagrange }}~~~~ /\!\!\!{\partial}/\!\!\!{\partial}\mathbf{A}^{\!\dagger} ~=~ /\!\!\!{\partial}\mathbf{F} ~=~0 \end{equation}\]
Thus the equations of motion, the Maxwell equations, are given by $/\!\!\!{\partial}\mathbf{F}=0$. If we work out the Lorentz invariant Lagrangian operator field we get.
\[\begin{equation} \mathcal{L} ~~=~~ \tfrac12\mathbf{F}\mathbf{F} ~~=~~ \tfrac12\left(E^2-B^2\right)\!I ~+~ \left(\vec{E}\cdot\vec{B}\right)i\gamma^5\end{equation}\]
The Lorentz scalar $\tfrac12(E^2-B^2)$ of the electromagnetic field is associated with the diagonal matrix $I$ and gives rise to the inhomogenious Maxwell equations. The pseudo scalar $\vec{E}\cdot\vec{B}$ of the electromagnetic field is associated with the pseudo scalar generator $i\gamma^5$ and gives rise to the homogenious Maxwell equations.
There is more information in the PDF here and there is a Mathematica file here.
This is all part of a much larger project: https://thephysicsquest.blogspot.com/