Solution (source code)

= Solution

Let $P^\mu=\partial\mathcal L/\partial(\partial_\mu h)$ and $Q^{\mu\nu}=\partial\mathcal L/\partial(\partial_\mu\partial_\nu h)$. Define the second coefficient with ordered repeated indices and $Q^{\mu\nu}=Q^{\nu\mu}$, so that
$$
\delta\mathcal L=\mathcal L_h\delta h+P^\mu\partial_\mu\delta h+Q^{\mu\nu}\partial_\mu\partial_\nu\delta h.
$$
This convention avoids double-counting a mixed derivative such as $h_{,12}=h_{,21}$. The first derivative term needs one <integration by parts>. The second derivative term needs two:
$$
\int_\Omega Q^{\mu\nu}\partial_\mu\partial_\nu\delta h
=\text{boundary terms}+\int_\Omega(\partial_\mu\partial_\nu Q^{\mu\nu})\delta h.
$$
Combining the terms gives
$$
\delta S=\int_\Omega\left[\mathcal L_h-\partial_\mu P^\mu+\partial_\mu\partial_\nu Q^{\mu\nu}\right]\delta h
+\int_{\partial\Omega}n_\mu\left[(P^\mu-\partial_\nu Q^{\mu\nu})\delta h+Q^{\mu\nu}\partial_\nu\delta h\right]d\Sigma.
$$
For compactly supported variations, or fixed $\delta h$ and its first derivatives at the boundary, the last integral vanishes. The <fundamental lemma of the calculus of variations> then yields the <Euler-Lagrange field equation with second derivatives>:
$$
\boxed{\frac{\partial\mathcal L}{\partial h}
-\partial_\mu\left[\frac{\partial\mathcal L}{\partial(\partial_\mu h)}\right]
+\partial_\mu\partial_\nu\left[\frac{\partial\mathcal L}{\partial(\partial_\mu\partial_\nu h)}\right]=0.}
$$
The plus sign of the final term comes from the two integrations by parts. In general this gives fourth-order equations, although the order in a particular coordinate depends on which second derivatives occur in the density.