I would like to confirm my answer since it is in direct contradiction from my book on physics. So suppose we have an action functional $$S(\phi) = \int_M \phi(\Delta + m^2)\phi dV,$$ where $M$ is a manifold and $\phi$ is in $C_c^\infty(U)$, the space of compactly supported smooth functions in an open set $U$.
So I get the variation $$S(\phi+t\psi) = S(\phi) + t(\int_M \psi(\Delta + m^2)\phi dV + \int_M \phi(\Delta + m^2)\psi dV) + t^2 S(\psi),$$ hence getting the functional derivative $$\int_M \psi(\Delta + m^2)\phi dV + \int_M \phi(\Delta + m^2)\psi dV.$$
But then the book claims that the Euler Lagrange equation is $(\Delta + m^2)\phi=0$. Is there some condition on $\psi$ that I don't know about or did I misunderstand some of the notations involved?
I wish to show that, under reasonable hypotheses,
$\displaystyle \int_m \psi(\nabla^2 + m^2) \phi \; dV = \int_m \phi(\nabla^2 + m^2) \psi \; dV. \tag 0$
First of all observe that
$\displaystyle \int_M \psi (\nabla^2 + m^2) \phi \; dV = \int_M \psi \nabla^2 \phi \; dV + \int_M m^2 \psi \phi \; dV, \tag 1$
and
$\displaystyle \int_M \phi (\nabla^2 + m^2) \psi \; dV = \int_M \phi \nabla^2 \psi \; dV + \int_M m^2 \phi \psi \; dV; \tag 2$
since evidently
$\displaystyle \int_M m^2 \psi \phi \; dV = \int_M m^2 \phi \psi \; dV, \tag 3$
we only need show that
$\displaystyle \int_M \psi \nabla^2 \phi \; dV = \displaystyle \int_M \phi \nabla^2 \psi \; dV; \tag 4$
we have the vector identity
$\nabla \cdot (fX) = \nabla f \cdot X + f \nabla \cdot X \tag 5$
for functions $f$ and vector fields $X$ (see here); if we set $f = \psi$ and $X = \nabla \phi$ this yields
$\nabla \cdot (\psi \nabla \phi) = \nabla \psi \cdot \nabla \phi + \psi \nabla \cdot \nabla \phi = \nabla \psi \cdot \nabla \phi + \psi \nabla^2 \phi, \tag 6$
or
$\psi \nabla^2 \phi = \nabla \cdot (\psi \nabla \phi) - \nabla \psi \cdot \nabla \phi, \tag 7$
which we may integrate over $M$:
$\displaystyle \int_M \psi \nabla^2 \phi \; dV = \int_M \nabla \cdot (\psi \nabla \phi) \; dV - \int_M \nabla \psi \cdot \nabla \phi \; dV. \tag 8$
We wish to apply the divergence theorem to the first integral on the right of (8); in order to so this, we must pause and carefully consider the regions in $M$ on which $\psi, \phi \ne 0$; that is, the support of of $\psi$ and $\phi$. We are given that $\phi$ (and I assume $\psi$) are in $C_c^\infty(U)$,
$\phi, \psi \in C_c^\infty(U), \tag 9$
$U \subset M$ open; it follows that the compact set
$K = \text{supp} \; \phi \cup \text{supp} \; \psi \subset U, \tag{10}$
and we require the existence of an open $\Omega \subset U$ with $\bar \Omega$ compact, $K \subset \Omega$, and
$\partial \Omega = \bar \Omega \setminus \Omega \tag{11}$
a smooth (at least $\mathcal C^1$), orientable, compact hypersurface in $M$. These conditions may seem somewhat arbitrary and or strange, but what I have tried to do in imposing them is simply ensure that we may perform, on the manifold $M$, the kind of calculations around the divergence theorem that we typically do in $\Bbb R^n$; e.g., assuming for the moment that in fact
$M = \Bbb R^n, \tag{12}$
then we see that we may take
$\Omega = B(R, 0) = \{ x \in \Bbb R^n \mid \vert x \vert < R \}, \tag{13}$
that is, $\Omega$ is the open ball of radius $R$ centered at the origin; then
$\bar \Omega = \bar B(R, 0) = \{ x \in \Bbb R^n \mid \vert x \vert \le R \}, \tag{14}$
and
$\partial \Omega = \bar \Omega \setminus \Omega = S(R, 0) = \{ x \in \Bbb R^n \mid \vert x \vert = R \} \tag{15}$
is the sphere of radius $R$ about the origin. Since $S(R, 0)$ is compact and has a well-defined outward pointing normal field $\mathbf n$, we may, for $R$ large enough, use the divergence theorem applied to the first integral on the right of (8):
$\displaystyle \int_{\Bbb R^n} \nabla \cdot (\psi \nabla \phi) \; dV = \int_{\bar B(R, 0)} \nabla \cdot (\psi \nabla \phi) \; dV = \int_{S(R, 0)} \psi \nabla \phi \cdot \mathbf n \; dA, \tag{16}$
where $dA$ is the volume measure on $S(R, 0)$; now by virtue of the hypothesis that both $\psi, \phi \in C_c^\infty(\Bbb R^n)$, we see that for $R$ sufficiently large,
$\displaystyle \int_{S(R, 0)} \psi \nabla \phi \cdot \mathbf n \; dA = 0, \tag{17}$
and so by (8) we have
$\displaystyle \int_{\Bbb R^n} \psi \nabla^2 \phi \; dV = \int_{\Bbb R^n} \nabla \cdot (\psi \nabla \phi) \; dV - \int_{\Bbb R^n} \nabla \psi \cdot \nabla \phi \; dV$ $= \displaystyle \int_{S(R, 0)} \nabla \cdot (\psi \nabla \phi) \; dV - \int_{\Bbb R^n} \nabla \psi \cdot \nabla \phi \; dV = - \int_{\Bbb R^n} \nabla \psi \cdot \nabla \phi \; dV. \tag{18}$
From this perspective, we see that they hypotheses we have placed upon $\Omega \subset M$ have been devised to allow a calculation similar to (18) to proceed on/within $M$; with such $\Omega$ as described above ca. equations (10)-(11) we may write
$\displaystyle \int_M \nabla \cdot (\psi \nabla \phi) \; dV = \int_{\bar \Omega} \nabla \cdot (\psi \nabla \phi) \; dV = \int_{\partial {\bar \Omega}} \psi \nabla \phi \cdot \mathbf n \; dA = 0, \tag{19}$
where again $\mathbf n$ is the outward-pointing normal on $\partial{\bar \Omega}$; we have used the fact that $\phi, \psi = 0$ on $\partial{\bar \Omega}$ in establishing this equation. In the light of (19) we see that (8) becomes
$\displaystyle \int_M \psi \nabla^2 \phi \; dV = -\int_M \nabla \psi \cdot \nabla \phi \; dV; \tag{20}$
we may reverse the roles of $\phi$, $\psi$ in (20) and arrive at
$\displaystyle \int_M \phi \nabla^2 \psi \; dV = -\int_M \nabla \phi \cdot \nabla \psi \; dV; \tag{21}$
by virtue of (20) and (21) in collusion we have
$\displaystyle \int_M \psi \nabla^2 \phi \; dV = \int_M \phi \nabla^2 \psi \, dV \tag{22}$
and thus, combining (3) and (22) we have
$S(\phi + t\psi) = S(\phi) + 2t \displaystyle \int_M \psi(\nabla^2 + m^2)\phi \; dV + t^2 S(\psi); \tag{23}$
we see that $S(\phi)$ is stationary precisely when
$\displaystyle \int_M \psi(\nabla^2 + m^2)\phi \; dV = 0, \; \forall \psi \in C_c^\infty(M), \tag{24}$
which shows the Euler-Lagrange equation must be
$(\nabla^2 + m^2)\phi = 0. \tag{25}$