This calculation is tedious.
You will need:
$$\delta R_{ab} = \nabla_{b}(\delta \Gamma^{c}_{ac}) - \nabla_{c}(\delta \Gamma^{c}_{ab})$$
$$\delta \Gamma^{c}_{ab} = \cfrac{1}{2}g^{cd}\Big(\nabla_{b}\delta g_{da} + \nabla_{a}\delta g_{db} - \nabla_{c} \delta g_{ab}\Big)$$
$$G_{ab} = R_{ab} - \cfrac{1}{2}g_{ab}R$$
First of all use the expression for the variation of the Christoffel symbols to obtain the variation of the Ricci tensor as covariant derivatives of variations of the metric tensor.
Then for the terms:
$$\delta R_{ab}\partial^{a}\phi \partial^{b}\phi \hspace{1.0cm} g_{ab}\delta(R)\partial^{a}\phi\partial^{b}\phi$$
you have to perform integration by parts to cancel total derivative terms. (see my answers here: Derivation of $f(R)$ field equations, problem with integration by parts , Metric field equations for the Jordan-Brans-Dicke action)
The calculation is not hard once you understand what you have to do, but it is going to take some time. I strongly encourage you to derive the field equations by hand, the satisfaction once you get the right answer does really worth the struggle.