From $$S=\int{d^4}x\sqrt{-g}f(R)$$ I want to deduce the $f(R)$ field equations which are: $$f'(R)R_{\mu\nu}-\frac{1}{2}f(R)g_{\mu\nu}-[\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\Box]f'(R)=\kappa T_{\mu\nu} ,$$
From the previous action I got so far $$\int{d^4x}\sqrt{-g}f'(R)R_{\mu\nu}\delta g^{\mu\nu}-\frac{1}{2}\sqrt{-g}g_{\mu\nu}f(R)\delta{g^{\mu\nu}}-\sqrt{-g}f'(R)[\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\Box]\delta{g^{\mu\nu}}.$$
It's clear that I got the first two terms of the field equations as the variation of the metric vanishes when we vary the action.
My problem is with the last two terms from the last equation.
I tried to integrate by parts twice and I got the following: \begin{eqnarray} f'(R)\nabla_\mu\nabla_\nu \delta g^{\mu\nu} &=&\nabla_\mu (f'(R)\nabla_\nu \delta g^{\mu\nu})-\nabla_\nu \delta g^{\mu\nu}\nabla_\mu f'(R) \\ &=&\nabla_\mu(\nabla_{\nu}\left(f'(R)\delta g^{\mu\nu})-\delta g^{\mu\nu}\nabla_{\nu}f'(R)\right)\\ & & - \nabla_{\nu}\left(\delta g^{\mu\nu}\nabla_{\mu}f'(R)\right)+\delta g^{\mu\nu}\nabla_\mu \nabla_\nu f'(R)\\ &=& \nabla_\mu \nabla_\nu(f'(R)\delta g^{\mu\nu})-2\nabla_{\mu}(\nabla_\nu f'(R)\delta g^{\mu\nu})+\delta g^{\mu\nu}\nabla _\mu \nabla_\nu f'(R) \end{eqnarray} And \begin{eqnarray} f'(R)g_{\mu\nu}\Box \delta g^{\mu\nu} &=& \nabla^{\alpha}(f'(R)g_{\mu\nu}\nabla_{\alpha}\delta g^{\mu\nu})-\nabla_{\alpha}\delta g^{\mu\nu} g_{\mu\nu}\nabla ^{\alpha}f'(R)\\ &=& \nabla^{\alpha}(\nabla_{\alpha}(f'(R)g_{\mu\nu}\delta g^{\mu\nu})-g_{\mu\nu}\delta g^{\mu\nu}\nabla_{\alpha}f'(R))\\& & -\nabla_{\alpha}(\delta g^{\mu\nu} g_{\mu\nu}\nabla^{\alpha}f'(R))+\delta g^{\mu\nu}g_{\mu\nu}\Box f'(R)\\ &=& \Box (f'(R)g_{\mu\nu}\delta g^{\mu\nu})-2\nabla_{\alpha}(g_{\mu\nu}\delta g^{\mu\nu}\nabla ^{\alpha}f'(R))+\delta g^{\mu\nu}g_{\mu\nu}\Box f'(R) \end{eqnarray} Now it's clear that the last terms are the wanted terms. The middle terms have $\delta g^{\mu\nu}$ in the boundary vanishes by Stokes's theorem. My problem is with the terms: $\nabla_\mu \nabla_\nu(f'(R)\delta g^{\mu\nu})$ and $\Box (f'(R)g_{\mu\nu}\delta g^{\mu\nu})$.
How can we resolve them to get the field quations.
p.s. Please correct me if I'm wrong in my calculations.
