I'm having several doubts about the procedure proposed by the Dirac-Bergmann algorithm in order to get the correct equations of motion of electrodynamics (Maxwell's equations).
Suppose I've already found the primary constraint, $\pi^0 \approx 0$, the secondary, $\partial_i \pi^i - j^0 \approx 0$, and verified that no more independent or inconsistent conditions arised, and checked that both constraints are first class.
After that, I'm confused about a lot of things. Perhaps it's better that I describe what I'm doing. So, I coupled both constraints to the hamiltonian density via multipliers:
\begin{equation} \mathcal{H}_{ext}=\frac{1}{4}F_{ij}F^{ij}-\frac{1}{2}\pi_i \pi^i -A_0\partial_i \pi^i + j^{\mu} A_{\mu} + \lambda_1 \pi^0 + \lambda_2 \partial_i \pi^i \end{equation}
where $F^{ij}$ is the usual electromagnetic tensor, $A_{\mu}=(V,\vec{A})$, $j^{\mu}=(\rho,\vec{j})$ and $\pi^\mu$ the canonical momentum conjugated to $A_\mu$. Then I proceeded to find the EOM:
\begin{equation} \dot{\pi^{\omega}}=\{\pi^{\omega},H_{ext}\}=-\frac{\delta H_{ext}}{\delta A_{\omega}}=\partial_i\frac{\partial \mathcal{H}}{\partial (\partial_i A_{\omega})}-\frac{\partial \mathcal{H}}{\partial A_{\omega}} \\ \end{equation} \begin{equation} \dot{A_{\omega}}=\{A_{\omega},H_{ext}\}=\frac{\delta H_{ext}}{\delta \pi^{\omega}}=\frac{\partial \mathcal{H}}{\partial \pi^{\omega}}-\partial_i\frac{\partial \mathcal{H}}{\partial (\partial_i \pi^{\omega})} \end{equation}
where the second equality comes trivially from the (evaluation of the) definition of the Poisson bracket. Making the calculations, I get the following
\begin{align} \dot{\pi^{0}}&=\partial_i F^{i0}+\partial_i \pi^i -j^0 \to \partial_i F^{i0} \approx 0 \\ \dot{\pi^k}&=\partial_i F^{ik}-j^k \end{align}
and
\begin{align} \dot{A_0}&=\lambda_1-\pi_0+\partial_0 A_0 - \partial_0 \lambda_2 \to \dot{A_0}\approx\lambda_1+\partial_0 A_0 - \partial_0 \lambda_2\\ \dot{A_i}&=\partial_i A_0 -\pi_i-\partial_i \lambda_2 \end{align}
So, my doubts:
(1) I assume my calculations are wrong because I'm not being capable of recognizing the inhomogeneous Maxwell's equations on my results as (I think) I should ($\partial_{\mu}F^{\mu \nu}=j^{\nu}$), and I can't find the error.
(2) Even though I know that this is part of the fundamentals of the theory, I can't figure out if I should use the hamiltonian with only the primary constraints (as my doubtful orientator says) coupled, or with all first class coupled or even with the complete extended hamiltonian (as I denoted here and in this case coincides with the case that has all first class constraints coupled), and, of course, don't know the reason why the right case is correct. I can't realize this from the textbooks that treat the subject on the matricial 'all-in-one' way.
I know that the right way to understand this remarks would be back up a little and start over (even from the prerequisite topics) in order to assimilate the theory until it's close from being intuitive, but I've been doing that exhaustively and getting only more confused.