You vary the action with brute force Variation for simplicity. I assume that you can handle the last two terms and the problem is the first one. Varying with respect to $\phi$ you'll obtain:
$$\Box \phi + \cfrac{1}{4}R\phi -m^2 \phi =0$$
Varying with respect to the metric you have to integrate by parts the first term.
$$\delta(\phi ^2 R) = \phi^2 \delta R = \phi^2 (g_{\mu\nu}\Box \delta g^{\mu\nu} - \nabla_{\mu} \nabla_{\nu}\delta g^{\mu\nu} + R_{\mu\nu}\delta g^{\mu\nu})$$
You want $\delta g^{\mu\nu} $ to be a multiplying factor so you integrate by parts twice and now the derivatives will act on $\phi^2$.
For the integration check my answer here: Derivation of $f(R)$ field equations, problem with integration by parts
Τhe final answer will be:
$$ \cfrac{1}{8}\left( g_{\mu\nu}\Box - \nabla_{\mu}\nabla_{\nu} + G_{\mu\nu}\right)\phi^2 - \cfrac{1}{2}\nabla_{\mu}\phi\nabla_{\nu}\phi + \cfrac{1}{4}g_{\mu\nu}\nabla^{\xi}\phi\nabla_{\xi}\phi + \cfrac{1}{4}g_{\mu\nu}m^2\phi^2=0$$
You will obtain $4$ Einstein's equations and one Klein-Gordon. Only two of them are independent though. You can of course plug in the metric and vary with respect to the scale factor. The result will be the same. If you want to vary the action by hand and then use a program to obtain the components of equations, varying with respect to the field surely contains less manipulations than varying with respect to the scale factor $a$.
