Proca Lagrangian manipulation How can I show that the Lagrangian density
$$\mathcal{L} = -\frac{1}{2}\partial_\alpha \varphi_\beta \partial^\alpha \varphi^\beta + \frac{1}{2} \partial_\alpha \varphi^\alpha \partial_\beta \varphi^\beta + \frac{\mu^2}{2}\varphi_\alpha \varphi^\alpha$$
for the real vector field $\varphi^\alpha$ leads to the field equations
$$[g_{\alpha\beta}(\square+ \mu^2)-\partial_\alpha\partial_\beta]\varphi^\beta=0$$
I have done through lot of work, but I was wondering if there are some tricks or efficient ways to do it? 
And second question, how can I show that the field satisfies the Lorenz condition $\partial_\alpha \varphi^\alpha = 0$?
Source: Mandl Shaw QFT problem 2.3
 A: To derive the field equations more quickly, consider that all the terms in the Lagrangian are "squares": a tensor contracted with itself on all indices. The variation of a square is $\delta(F_{ab} F^{ab}) = 2 (\delta F_{ab}) F^{ab}$ in analogy with $d/dx\; f^2 = 2ff'$. Using this, we have for the variation of the Lagrangian $$\delta \mathcal L = - (\delta \partial_a \varphi_b) \partial^a \varphi^b + (\delta \partial_a \varphi^a) (\partial_b \varphi^b) + \mu^2 (\delta \varphi_a) \varphi^a.$$
We can swap variation and partial derivatives and find that the variation in the action is $$\delta S =\int \delta \mathcal L = \int (\delta \varphi_a )\partial_\nu \partial^\nu \varphi^a - (\delta \varphi^a)(\partial_a \partial_b \varphi
^b) + \mu^2 (\delta \varphi_a) \varphi^a$$
where we used integration by parts on the first two terms and relabeled the summation indices in the first. We see that the variation is $$\delta S = \int (\delta \varphi^a) (g_{ab}\square  \varphi^b - \partial_a \partial_b + \mu^2 g_{ab} )\varphi^b$$
which vanishes for an arbitrary variation $\delta \varphi^a$ if and only if the field equation $$(g_{ab}(\square + \mu^2) - \partial_a\partial_b ) \varphi^b = 0$$
holds. For the second part: apply $\partial^a$ to the equation. This yields $$(\partial_b (\square + \mu^2) - \partial^a \partial_a \partial_b) \varphi^b = \mu^2 \partial_b \varphi^b + (\partial_b\square - \square \partial_b) \varphi^b = 0,$$
which implies the Lorentz condition $\partial_b\varphi^b=0$ since the term with brackets vanishes by commutativity of partial derivatives.
