Hamiltonian field equations constraints 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?
 A: OP wrote (v3):

Is there anything in particular I should be careful of?

Yes. Watch out for secondary constraints, cf. e.g. this Phys.SE post. 
Below follows a brief partial derivation. 


*

*Let Greek letters $\mu,\nu,\ldots$ denote spacetime indices, while Roman letters $i,j,\ldots$ denote only spatial indices. The Lagrangian density
$$ {\cal L}~=~ -\frac{1}{2}\partial_{\mu}\phi^{\nu}\partial^{\mu} \phi_{\nu} +\frac{1}{2}(\phi^{\mu}_{,\mu})^2+\frac{1}{2}m^2\phi^{\mu}\phi_{\mu}$$
$$~=~\frac{1}{2}\dot{\phi}^i\dot{\phi}^i+ \dot{\phi}^0\phi^i_{,i} +\frac{1}{2}(\phi^i_{,i})^2-{\cal V},\tag{A}$$
with potential density
$$-{\cal V}~:=~\frac{1}{2}\partial_i\phi^{\mu}\partial_i\phi_{\mu} +\frac{1}{2}m^2\phi^{\mu}\phi_{\mu}, \tag{B}$$
and with Minkowski metric $\eta_{\mu\nu}={\rm diag}(+1,-1,-1,-1)$. 

*The momentum is equal to the velocity
$$ \pi_i~=~\dot{\phi}^i,\tag{C}$$
except for $\pi_0$. 

*The primary constraint reads 
$$ \chi^1~:=~\pi_0- \phi^i_{,i}~\approx~0.\tag{D}$$

*The Hamiltonian density becomes
$$ {\cal H}~=~-\frac{1}{2}\pi_{\mu}\pi^{\mu}+{\cal V}. \tag{E}$$

*The CCRs read
$$ \{\phi^{\mu}({\bf x}),\pi_{\nu}({\bf y})\}_{PB}~=~\delta^{\mu}_{\nu}~\delta^3({\bf x}-{\bf y}),\tag{F}$$
and all other vanish.

*A secondary constraint 
$$ \chi^2~:=~ m^2\phi^0 - \phi^0_{,ii} - \pi_{i,i}~\approx~0\tag{G}$$
arises because $\{\chi^1,H\}_{PB}$ is not proportional to $\chi^1$. Without a secondary constraint, the time evolution would violate the primary constraint (D). One can check that there doesn't appear a tertiary constraint. 

*The Poisson bracket between the two constraints reads
$$ \{\chi^2({\bf x}),\chi^1({\bf y})\}_{PB}
~=~m^2\delta^3({\bf x}-{\bf y}).\tag{H} $$ 
In other words, the two constraints are first (second) class if $m^2=0$ ($m^2\neq 0$), respectively. 

*First class constraints generate gauge symmetry. 

*For second class constraints, the Dirac bracket 
$$\{\phi^0({\bf x}), \pi_0({\bf y})\}_{DB}~=~\frac{\partial_i\partial_i}{m^2}\delta^3({\bf x}-{\bf y}), $$
$$\{\phi^i({\bf x}), \pi_j({\bf y})\}_{DB}~=~\left(\delta^i_j-\frac{\partial_i\partial_j}{m^2}\right)\delta^3({\bf x}-{\bf y}), $$
$$\{\phi^0({\bf x}), \phi^i({\bf y})\}_{DB}~=~-\frac{\partial_i}{m^2}\delta^3({\bf x}-{\bf y}),\qquad \text{etc},\tag{I}$$
is necessary.
