Generally covariant Klein-Gordon equation Consider a 4-dimensional smooth manifold $M$ on which there is a Lorentzian metric $g_{ab}$ and a function $\phi$ satisfying the following two equations (in abstract index notation):
\begin{equation}
g^{ab}\nabla_{a}\nabla_{b}\phi - m^{2}\phi = 0,
\qquad
R_{abcd} = 0.
\end{equation}
I've been told that the covariance group of these equations is $\mathrm{Diff}(M)$, but I am having a hard time seeing it for myself.  From what I understand, this means that for any $\psi \in \mathrm{Diff}(M)$, $\psi^{*}\phi$ (and $\psi^{*}g_{ab}$?) is a solution of these equations just if $\phi$ (and $g_{ab}$?) is.  It's straightforward that $\psi^{*}R_{abcd} = 0$, but I am having a hard time showing that $\psi^{*}(g^{ab}\nabla_{a}\nabla_{b}\phi) = g^{ab}\nabla_{a}\nabla_{b}(\psi^{*}\phi)$.  Is this true?  If so, how would I go about showing that?
 A: The answer given by @YuriVyatkin is the right answer to my question, but I thought I'd write out more fully how it answers my question.
Usual statements of general covariance are something like "the form of the equations is the same in every frame" or "the diffeomorphism group is the covariance group of the theory."  While I'm sure these are true, they are too imprecise for me to understand (or involve terms like "form" that I can't find the mathematical definition of), which is why I wanted some mathematical statement that corresponds to what they mean.  Sometimes physicists will say that if you have a tensor equation, then you automatically have general covariance.  I could accept that as a definition of general covariance, but it wasn't clear to me how that's related to the other two characterizations.  I've settled on something that I find satisfying, which boils down to the naturality of the pullback.
The basic statement of general covariance that seems right to me is from Anderson's Principles of Relativity Physics, which in a slightly modernized form is something like:
A group $G$ is a covariance group of a theory if


*

*There is a faithful representation $\rho$ of $G$ on the state space of the theory, and

*If $\langle O_{1}, O_{2}, \dots \rangle$ is a solution of the equations of the theory (where $O_{i}$ are some quantities), then for all $g \in G$, $\langle \rho(g)O_{1}, \rho(g)O_{2}, \dots \rangle$ is a solution as well.
So consider a massive scalar Klein-Gordon field on flat spacetime:
$$
g^{ab}\nabla_{a}\nabla_{b}\phi - m^{2}\phi = 0, \qquad
R_{abcd} = 0.
$$
The state space of the theory is pairs $\langle g_{ab}, \phi \rangle$ of a
Lorentzian metric and a scalar field, and there is a faithful representation of $\mathrm{Diff}(M)$ on these pairs given by
$$
\rho(\psi)\colon \langle g_{ab}, \phi \rangle \mapsto \langle (\psi^{-1})^{*}g_{ab}, (\psi^{-1})^{*}\phi \rangle
$$
for each $\psi \in \mathrm{Diff}(M)$.  We're allowed to pull back $g_{ab}$
    because we're only concerned with diffeos.  (And actually, we can get rid
    of the inverses by noting that what we're really doing is pushing both
    things forward.  But I'll stick with talking in terms of pullbacks.)  The
    result of applying $\rho(\psi)$ is a pair of a Lorentzian metric and a
    scalar field, so we're still in the state space of the theory, and since
    the pushforward is linear, it really is a representation.  Faithfulness is
    also pretty direct.  So now we just need to show that $\langle \psi^{*}g_{ab}, \psi^{*}\phi \rangle$ is a solution to the equations for all $\psi^{-1} \in \mathrm{Diff}(M)$.  
This is what was really hanging me up, since I wasn't sure how to see this.  I now realize that part of my problem was the fact that $\nabla_{a}$ was involved, and I didn't know how pullbacks interacted with it.  So consider a simpler theory for a moment, of the form
$$
G_{ab} = 8\pi T_{ab}
$$
To show that $\langle \psi^{*}G_{ab}, \psi^{*}T_{ab} \rangle$ is a solution is straightforward:
$$
(\psi^{*}G)_{ab} - 8\pi(\psi^{*}T)_{ab}
=
\psi^{*}(G_{ab} - 8\pi T_{ab})
=
\psi^{*}0
=
0
$$
which follows directly from the linearity and naturality of the pullback with respect to all the tensor operations (scalar multiplication, tensor addition, etc.).  This example makes it clear why tensor equations are key: since the pullback is natural, you'll always be able to bring it all outside the equation, and the right side will always vanish.  So if you have a tensor equation, then you get diffeomorphism covariance for free.
Okay, back to the Klein-Gordon equation.  The tensor equations involve 5 doodads ($g_{ab}$, $\nabla_{a}$, $\phi$, $m$, and $R_{abcd}$), but only two show up in the state pairs.  Which is fine, since $\nabla_{a}$ is assumed to be the Levi-Civita covariant derivative, so it and $R_{abcd}$ are determined by $g_{ab}$.  And since $m^{2}$ is just a constant function, it's invariant under the pullback.  Since the pullback of the Levi-Civita connection is the Levi-Civita connection, and the pullback of a flat metric is a flat metric, the second equation is automatically satisfied, and for the first we have
$$
(\psi^{*}g)^{ab}(\psi^{*}\nabla)_{a}(\psi^{*}\nabla)_{b}(\psi^{*}\phi)
-
(\psi^{*}m^{2})(\psi^{*}\phi)
=
\psi^{*}(g^{ab}\nabla_{a}\nabla_{b}\phi - m^{2}\phi)
=
\psi^{*}0
=
0.
$$
So there you go, $\langle \psi^{*}g_{ab}, \psi^{*}\phi\rangle$ is really a solution.  I think the biggest problem I was having was seeing what the importance of the pullback was.  It's the reason that tensor equations give us general covariance, but I've never seen someone make this connection explicit.  I think it's also what people mean when they talk about the "form" of the equations being preserved under coordinate transformations.  I still am not totally sure what "form" means here, but it's not really important, since I can just translate that to naturality when I hear it.
Also, just to quickly answer the last question from my post:  it's not true in
general that 
$
g^{ab}\nabla_{a}\nabla_{b}(\psi^{*}\phi) =
g^{ab}\nabla_{a}\nabla_{b}\phi
$.
Two ways to see this:  the way that I first realized it is that
$\nabla_{a}\nabla_{b}\phi = {\star}\, d\star d\phi$, and the pullback only
commutes (up to a sign) with the Hodge star when $\psi$ is an isometry.  Since
there are two stars, any pullback by an isometry will satisfy
$
g^{ab}\nabla_{a}\nabla_{b}(\psi^{*}\phi) =
g^{ab}\nabla_{a}\nabla_{b}\phi
$.
Alternatively, note that since $\nabla_{a}$ is determined by $g_{ab}$, this
equation is just the demand that $(\psi^{*}g)_{ab} = g_{ab}$; i.e., that $\psi$
is an isometry.
A: Wouldn't it better to have $\psi^{*}(g^{ab}\nabla_{a}\nabla_{b}\phi) = (\psi^{*}g)^{ab}(\psi^{*}\nabla)_{a}(\psi^{*}\nabla)_{b}(\psi^{*}\phi)$?
Then this is just the naturality of the metric and the connection.
