Change of variables in integration on manifolds: what's wrong? The problem is this: Let $f:M^m\to\mathbb{R}$ a continuous function with compact support and $h:N\to M$ a diffeomorphism between Riemannian compact manifolds. Then
$$\int_{M}fdV^M=\int_N f\circ h|\det dh|dV^N$$
The proof using differential forms is easy, but I'm trying to prove this without using forms, just the definition.
"Proof":
Let $\Omega$ be the support of $f$, $\Theta=h^{-1}(\Omega)$ the support of $f\circ h$, $\textbf{x}_\alpha:U_\alpha\to M$ and $\textbf{y}_\alpha:V_\alpha\to N$ parametrizations and $\{\phi_\alpha\}_{\alpha\in A}$ a partition of unity subordinate to $\{U_\alpha\}_{\alpha\in A}$. Let $g^M$ the metric of $M$ and $g^N$ the metric of $N$.Then
$$\int_M fdV^M=\sum_{\alpha}\int_{\textbf{x}_\alpha^{-1}(\Omega)}\phi_\alpha f\circ\textbf{x}_\alpha\sqrt{\det g^M}dx^M$$
Denoting $\Omega_\alpha=\textbf{x}_\alpha^{-1}(\Omega)\subset U_\alpha$,$\Theta_\alpha=\textbf{y}^{-1}_\alpha(\Theta)\subset V_\alpha$ and $T_\alpha=\textbf{x}_\alpha^{-1}\circ h\circ\textbf{y}_\alpha:V_\alpha\to U_\alpha$ we will have
$$\int_M fdV^M=\sum_{\alpha}\int_{\Omega_\alpha}\phi_\alpha f\circ\textbf{x}_\alpha\sqrt{\det g^M}dx^M=\sum_{\alpha}\int_{T_\alpha(\Theta_\alpha)}\phi_\alpha f\circ\textbf{x}_\alpha\sqrt{\det g^M}dx^M$$
by the variable change theorem,
$$\int_M fdV^M=\sum_{\alpha}\int_{\Theta_\alpha}(\phi_\alpha\circ T_\alpha)(f\circ\textbf{x}_\alpha\circ T_\alpha)\sqrt{\det g^M\circ T_\alpha}|\det dT_\alpha|dy^N$$
How
\begin{eqnarray*}
\sqrt{\det g^M\circ T_\alpha}&=&\mathrm{vol}\ P\left[\frac{\partial}{\partial x_1},\cdots,\frac{\partial}{\partial x_m}\right] \\
&=&\mathrm{vol}\ P\left[dh\left(\frac{\partial}{\partial y_1}\right),\cdots,dh\left(\frac{\partial}{\partial y_m}\right)\right] \\
&=&\mathrm{vol}\ P\left[\frac{\partial}{\partial y_1},\cdots,\frac{\partial}{\partial y_m}\right]|\det dh| \\
&=&\sqrt{g^N}|\det dh|
\end{eqnarray*}
and
$$\det dT_\alpha=\det dh$$
We have
\begin{eqnarray*}
\int_M fdV^M&=&\sum_{\alpha}\int_{\Theta_\alpha}(\phi_\alpha\circ T_\alpha)(f\circ\textbf{x}_\alpha\circ T_\alpha)\sqrt{\det g^N}|\det dh|^2dy^N \\
&=&\sum_{\alpha}\int_{\textbf{y}_\alpha^{-1}(\Theta)}(\phi_\alpha\circ T_\alpha)(f\circ h\circ \textbf{y}_\alpha)\sqrt{\det g^N}|\det dh|^2dy^N \\
&=&\int_N f\circ h|\det dh|^2dV^N
\end{eqnarray*}
Question: where is the error? It must probably be something stupid but I can't see it.
 A: Suppose $(M,g)$ is a Riemannian manifold, and $h:N\to M$ a diffeomorphism. Then, we can consider the Riemannian manifold $(N,h^*g)$. Then, for any smooth function $f:M\to\Bbb{R}$, we can write the change of variables formula as
\begin{align}
\int_Mf\,dV_g &=\int_{N}f\circ h\, dV_{h^*g}\tag{i}
\end{align}
This should remind you very much of the differential form version where $\int_M\omega = \int_{h(N)}\omega = \int_Nh^*\omega$. Here, you can think of $dV_g$ as the Riemannian volume form on $M$ induced by the metric $g$ (provided $M$ is orientable, in which case the choice of volume form provides the orientation). Or you can also think of this as a measure (which you can do even when $M$ is not orientable).
As I mentioned in the comments, writing $\det dh$ abstractly as you've done is incorrect (at the very least very bad practice) because (pointwise) $dh$ is a linear transformations between different vector spaces. Now, you seem to be starting with two Riemannian manifolds $(M,g^M)$ and $(N,g^N)$, and then considering a diffeomorphism $h:N\to M$. In this case, one can also write down a change of variables formula:
\begin{align}
\int_Mf\,dV_{g^M} &=\int_{N}f\circ h\cdot \left(\dfrac{dV_{h^*(g^M)}}{dV_{g^N}}\right)\, dV_{g^N}\tag{ii}
\end{align}
Here, the bracketed term should be interpreted as a Radon-Nikodym derivative of the two measures $dV_{h^*(g^M)}$ and $dV_{g^N}$ (which one can show are positive $\sigma$-finite measures which are mutually absolutely continuous, hence the Radon-Nikodyn theorem can be applied). This is the correct factor which should appear (very roughly speaking, it's like writing $dy=\dfrac{dy}{dx}dx$ inside the integral sign, or writing $d^ny=\left|\det\frac{\partial y}{\partial x}\right|\,d^nx$ in the multivariable setting). As you can see from the notation itself, this Radon-Nikodym derivative involves all three things: $h,g^M$ and $g^N$, and this takes into account the "change of measure" factor (which you incorrectly just wrote as $\det dh$). Of course, if one works things out in local coordinates, one will get some sort of determinant of a coordinate representation of derivatives of $h$, but in addition, one will also get extra factors involving the determinants of $g^M, g^N$ in appropriate coordinates.

As you can see (ii) is not the nicest way of writing things down; (i) is the nice way of doing things. Now, let us forget about (ii) and try to prove (i), atleast locally. So, let's forget the excess baggage of a partition of unity and work with a single chart (because like I said, once we do this, we just put in a few extra summations and we're done).
To this end, let $(U,x)$ be a chart on $M$, and $(V,y)$ a chart on $N$ such that $h(V)\subset U$ (for my own sanity, I'm working with charts... i.e the function $x:U\to x[U]\subset\Bbb{R}^n$ goes from the manifold to the Euclidean space, rather than the other way around). Also, I shall denote the metric on $M$ simply by $g$, and I consider the pull-back metric $h^*g$ on $N$. We now want to show that
\begin{align}
\int_{h(V)}f\,dV_g =\int_{V}f\circ h\,dV_{h^*g}.
\end{align}
Ok so let's start on the right. We have
\begin{align}
\int_{V}f\circ h\,dV_{h^*g} &:= \int_{y(V)}(f\circ h)\circ y^{-1}\cdot
\sqrt{\left| \det [(h^*g)_{ij}]\right|\circ y^{-1}}\, dy\tag{$*$}
\end{align}
Let's calculate what the entries of the matrix are at a point $b\in y(V)$:
\begin{align}
(h^*g)_{ij}(y^{-1}(b))&:=(h^*g)\left(
\frac{\partial}{\partial y^i}\bigg|_{y^{-1}(b)},
\frac{\partial}{\partial y^j}\bigg|_{y^{-1}(b)}\right)\\
&=g\left[h_*\left(\frac{\partial}{\partial y^i}\bigg|_{y^{-1}(b)}\right),
h_*\left(\frac{\partial}{\partial y^j}\bigg|_{y^{-1}(b)}\right) \right]\\
&= \frac{\partial (x^{\mu}\circ h)}{\partial y^i}\bigg|_{(x\circ h\circ y^{-1})(b)}\cdot 
\frac{\partial (x^{\nu}\circ h)}{\partial y^j}\bigg|_{(x\circ h\circ y^{-1})(b)}\cdot
g_{\mu\nu}(x\circ h\circ y^{-1}(b))
\end{align}
The efficient way of writing this down is that the matrices are related via a congruence (I'm omitting the points of evaluation)
\begin{align}
[h^*g] &= [D(x\circ h\circ y^{-1})]^t\cdot [g_{\mu\nu}]\cdot [D(x\circ h\circ y^{-1})]
\end{align}
(This should hopefully be familiar from linear algebra: the matrix representations of a bilinear map are related by congruence). Now, since a matrix and its transpose have the same determinant, it follows that the determinants are related as:
\begin{align}
\sqrt{\left|\det (h^*g)_{ij}(y^{-1}(b))\right|}&=
\sqrt{\left|\det g_{\mu\nu}((x\circ h\circ y^{-1})(b))\right|}\cdot
\left|\det D(x\circ h\circ y^{-1})_b\right|
\end{align}
We now have everything we need to use the standard change of variables in the RHS of $(*)$:
\begin{align}
\text{RHS of $(*)$} &=
\int_{y(V)} (f\circ x^{-1})\circ (x\circ h\circ y^{-1})\cdot
\sqrt{\left|\det g_{\mu\nu}\right|}\circ (x\circ h\circ y^{-1})\cdot
\left|\det D(x\circ h\circ y^{-1})\right|\,dy\\
&=\int_{x(h(V))}f\circ x^{-1}\cdot \sqrt{\left|\det g_{\mu\nu}\right|}\,dx\\
&:=\int_{h(V)}f\, dV_g.
\end{align}
This completes the proof.
