Hamilton's equations of motion on Dirac's formalism 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.
 A: It is, a priori, completely correct to add both primary and secondary constraints to the Hamiltonian density by Lagrange multipliers. What is not correct is how you determined the equations of motion:


*

*There is no "$F^{i0}$" in the Hamiltonian theory! It is called $\pi^i$ there and it is not dependent on $\partial_i A_0$, it is an independent canonical variable! I don't know how you got the $F^{i0}$ in your expression for $\dot{\pi}^0$, but it's not there. The correct equation of motion is just the secondary constraint
$$ \dot{\pi}^0 = \partial_i \pi^i - j^0\tag{1a}$$

*The equation of motion
$$ \dot{\pi}^k = \partial_i F^{ik}\tag{1b} - j^k$$
is correct.

*The equation of motion for $A_0$ is just
$$ \dot{A}_0 = \lambda_1\tag{1c}$$
and I don't know where your other terms come from. The only term in the Hamiltonian that depends on $\pi^0$ is $\lambda_1\pi^0$, and no term depends on $\partial_i \pi^0$.

*The equation of motion for $A_i$,
$$ \dot{A}_i = -\pi_i + \partial_i A_0 - \partial_i \lambda_2 \tag{1d}$$
is correct.


Maxwell's equations are now obtained as follows: 


*

*Identifying $E^i = -\pi^i$, we see that $\nabla\cdot E = \rho$ is not an equation of motion for $A$ or $\pi$, but a constraint or the equation of motion for $\lambda_2$, whichever you prefer. In particular, eq. (1a) has no dynamical content.

*Identifying $B_i \propto \epsilon_{ijk}F^{jk}$, eq. (1b) is $\nabla\times B - \dot{E} = j$.

*Eq. (1c) has no dynamical content, as it involves only Lagrange multipliers. We may make the gauge choice $\lambda_1 = 0$ to render $A_0$ constant in time.

*Eq. (1d) becomes $\pi_i = \partial_i A_0 - \dot{A}_i$ upon the gauge choice $\lambda_2 = 0$, which is just the equation $E = \nabla\phi -\partial_t A$.


If you are wondering where the other half of Maxwell's equations went - it's in the definition of the $F$! Both $\nabla\cdot B = 0$ and $\nabla\times E + \partial_t B = 0$ are direct consequences of defining $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$, since it follows from this that $\partial_\sigma \epsilon^{\sigma\rho\mu\nu}F_{\mu\nu} = 0$. $\sigma = 0$ gives $\nabla\cdot B = 0$ and $\sigma = i$ gives the $i$-th component of $\nabla\times E + \partial_t B = 0$.
