Proof of Liouville`s Theorem
From TorontoMathWiki
Consider the autonomous ODE $$ y'(t)=f(y(t)) $$ where $f \in C^{\infty}(\mathbb{R}^n, \mathbb{R}^n)$. The flow $F(t,y_0)$ associated with this equation is defined as follows. Let $y(t)$ denote the (unique) solution to this equation with initial condition $y(0) = y_0$ and $t \in I_x$, where $I_x$ is the maximal possible interval about $0$ where the solution exists. Then $F(t,y_0) : = y(t)$ for $t \in I_x$. Notice that $x \mapsto F(t,x)$ is smooth and injective (Why?). For a bounded open set $U \subset \mathbb{R}^n$, set $$V(t) = \text{Vol}\,(F(t, U)) : = \int_{F(t,U)}1 \,dx$$
Theorem: If $U\subset \mathbb{R}^n$ is a bounded open set, then $$\frac{d\, V}{dt} = \int_{F(t,U)} \operatorname{div} f(x) \, dx$$
Proof:
Notice that $F(t_0 + h, U) = F(h, F(t_0,U))$ provided that the right hand side is defined, which will be the case if $|h|$ is sufficiently small.
By a change of variables

We claim that for $x \in F(t_0, U)$, $$\frac{\partial F}{\partial x}(t,x) = I + tDf(x) + o(t)$$
if $|t|$ is sufficiently small. Indeed, note that $$F(t,x) = x + \int_{0}^t f(F(s,x))ds,$$ and after taking partial derivatives with respect to $x$ we get $$\frac{\partial F}{\partial x}(t,x) = I + \int_{0}^t \frac{\partial}{\partial x}f(F(s,x)) ds,$$ so that $$\frac{\partial}{\partial t} \frac{\partial F}{\partial x}(t,x) = \frac{\partial}{\partial x}f(F(t,x))$$ Evaluating at $t=0$ and using the chain rule we get $${\frac{\partial}{\partial t} \frac{\partial F}{\partial x}(t,x)}\Big|_{t=0} = Df(F(0,x))\cdot \frac{\partial F}{\partial x}(0,x) = Df(x)$$ The claim follows from this and from the definition of the derivative.
From the claim we deduce that $\operatorname{det} \left(\frac{\partial F}{\partial x}(t,x) \right) = \operatorname{det} (I + t\,Df(x)) + o(t)$. In the expansion of the determinant, the only summand that does not contribute an $o(t)$ term is $$ \prod_{i=1}^n(1 + t \frac{\partial f_{i}}{\partial x_{i}})= 1 + t \sum_{i=1}^n \frac{\partial f_{i}}{\partial x_{i}} + o(t)= 1 + t \operatorname{div} f(x) + o(t) $$
It follows that $$V(t_0+h) = \int_{F(t_0,U)} 1 + h\operatorname{div} f(x) + o(h) dx = V(t_0) + h \int_{F(t_0,U)} \operatorname{div} f(x) \,dx + o(h)$$ and this establishes the theorem from the definition of the derivative.
Corollary: The flow of a Hamiltonian system preserves volume since the corresponding $f = (-H_q, H_p)$ satisfies $\operatorname{div} f \equiv 0$.
Reference: Teschl, Gerald: Ordinary Differential Equations and Dynamical Systems, 2004. [1]