Let's consider the Lagrangian
$$\mathcal{L}~=~-\frac{1}{2}(\partial_\mu\phi^\nu)^2+\frac{1}{2}(\partial_\mu\phi^\mu)^2+\frac{1}{2}m^2\phi_\mu \phi^\mu,$$
with Minkowski metric $\eta_{\mu\nu}={\rm diag}(+1,-1,-1,-1)$
The equations of motion are then
$$-\partial_\nu\partial^\nu\phi_\mu + \partial_\mu\partial^\nu\phi_\nu-m^2\phi_\mu~=~0.$$
Taking the four divergence we find $\partial_\mu\phi^\mu=0$ so we can reduce the equations of motion to $$(\partial_\nu\partial^\nu +m^2)\phi_\mu~=~0.$$
This is all fine. However when I convert to the Hamiltonian formalism things don't quite work out. The conjugate momenta is given by
$$\pi_\mu~=~\frac{\partial\mathcal{L}}{\partial(\partial_0\phi^\mu)}~=~-\partial^0\phi_\mu + \partial^\nu\phi_\nu\delta^0_\mu$$
from which we get the constraint relation
$$\pi_0-\partial^i\phi_i~=~0$$
with $i=1,2,3$.
The Hamiltonian $\mathcal{H}$ is given by
$$\mathcal{H}~=~\pi_\mu\partial_0\phi^\mu-\mathcal{L}$$
and the Hamiltonian equations are
$$\frac{\partial \mathcal{H}}{\partial\pi_\mu}~=~\partial_0\phi^\mu$$ and $$\frac{\partial \mathcal{H}}{\partial\phi_\mu}~=~-\partial_0\pi^\mu.$$
Now when I explicitly write out the Hamiltonian and include the constraints etc the expression is getting quite complicated. Then calculating the eom is producing something not very nice which doesn't match my original equation of motion. EDIT: Just realised it shouldn't anyway since its one second order PDE vs two first order PDEs. The eom still aren't coming out as they should though.
Is there anything in particular I should be careful of?