First let me state the following elementary auxiliary Lemma from linear algebra
Lemma:
Let $\Delta:\mathbb{R}^{n^2}\longrightarrow\mathbb{R}$ be the
determinant function, i.e.
$$\Delta(\alpha_{11},\ldots,\alpha_{n1},\ldots,\alpha_{1n},
\ldots,\alpha_{nn})^{\top}
= \det[(\alpha_{ij})]$$
where $(\alpha_{ij})$ is the $n\times n$--matrix whose $ij$--th component is
$\alpha_{ij}$. Then,
$$\Delta_\alpha=
\frac{\partial \Delta}{\partial\alpha}=
(W_{11}\ldots,W_{n1},\ldots,W_{1n},\ldots,W_{nn})$$
where $W_{ij}$ is the $ij$--th cofactor of the matrix $(\alpha_{ij})$.
The proof of this Lemma an simple exercise about computing determinants using the cofactor formula.
Here is a proof of the statement of the OP:
Given the solution $\phi(t;{\bf x})=(\phi^1(t;{\bf x}),\ldots,\phi^n(t;{\bf x}))^\top$ to the initial value problem
$$ \dot{\mathbf{y}}(t)=f(t,\mathbf{y}(t)),\qquad \mathbf{y}(0)=\mathbf{x}$$
we use the notation
$\phi^{i}_{x_j}(t;\mathbf{x})= \frac{\partial\phi^i}{\partial x_j}(t;\mathbf{x})$.
Using the auxiliary Lemma above along with the chain rule we obtain
$$
\begin{align}
\dot{W}&= \sum_i W_{i1}\dot{\phi}^{i}_{x_1} +\cdots+ \sum_i
W_{in}\dot{\phi}^{i}_{x_n}\\
&=\sum_{ij} W_{ij}\dot{\phi}^{i}_{x_j}
\tag{1}\label{chain}
\end{align}
$$
where $W_{ij}$ is the $ij$--th cofactor of the matrix
$\left(\phi^i_{x_j}\right)$. It is easy to cheach that
$\phi_{\bf x}(t;{\bf x})=\frac{\partial\phi}{\partial
{\bf x}}(t;{\bf x})$ satisfies the variational equation
$$\begin{align}
\begin{matrix}
\dot{\phi}_{\bf x}(t;{\bf
x})&=&f_{\bf{x}}(t,\phi(t;{\bf{x}}))\phi_{\bf x}(t;{\bf x})\\
\phi_{\bf x}(0;{\bf x})&=&I
\end{matrix}
\tag{2}\label{vareq}
\end{align}
$$
Substituting \eqref{vareq} in \eqref{chain} and recalling the fact that
the determinant of a matrix that has two identical columns is zero, we obtain
$$
\begin{align}
\dot{W}(t)&=\sum_{ijk} W_{ij}(t)
f^i_{x_k}(t,\phi(t;{\bf x}))\phi^k_{x_j}(t;{\bf x})\\
&= \sum_{ki} \left(f^i_{x_k}(t,\phi(t;{\bf x})\right)
\sum_j W_{ij}(t)\phi^k_{x_j}(t;{\bf x})\\
&=\sum_i f^i_{x_i}(t,\phi(t;{\bf x})
\sum_j W_{ij}(t)\phi^i_{x_j}(t;{\bf x})\\
&= \sum_i f^i_{x_i}(t,\phi(t;{\bf x})) W(t) =
W(t)\,\left(\nabla_{\bf x}\cdot
f\right)(t,\phi(t;{\bf x}))
\end{align}
$$
Comments:
- The statement of the OP has the Wronskian equation for linear as a particular case. Indeed, for a non--autonomous linear system
$$
\dot{x}=A(t) x\quad x(0)=x_0,
$$
where $A(t)\in L(\mathbb{R}^n,\mathbb{R}^n)$ is continuously
differentiable in $t$, $f(t,\mathbf{x})=A(t)\mathbf{x}$ and so, $(\nabla_x\cdot f)(t,\mathbf{x})=\operatorname{Tr}(A(t))$ where $\operatorname{Tr}(A(t))= \sum_j a_{jj}(t)$, the trace of $A(t)$. Then $W$ satisfies
$$
\dot{W}=\operatorname{Tr}(A(t))\, W
$$ and so,
$$
W(t)=W(0) \exp\Big(\int^t_0\operatorname{Tr}(A(s))\,ds\Big)
$$
- The statement can be used to prove Liouville's theorem. The uniqueness theorem of solutions to initial value problem
$$ \begin{align}
\dot{y}=f(t, y),\quad y(s)=x\tag{3}\label{three}
\end{align}
$$
shows that the solutions to $\eqref{three}$
satisfy the following flow property:
If $\phi_{t,s}(y)=\phi(t;s,y)$ denotes the solution
in $t$ initial conditions $\phi(s;s,y)=y$, then
a. $\phi_{t,s}\circ\phi_{s,r}(y)=\phi_{t,r}(y)$ for all
$r,\,s,\,t\in I$ and $y\in \Omega$.
b. $\phi_{t,t}(y)=y$ for all $t\in I$ and $y\in\Omega$.
For simplicity, we assume that any solution
starting in $\Omega$ exists for all times, that is $I=\mathbb{R}$.
Suppose that $D(0)\subset\Omega$ has a finite volume $v(0)$ in $\mathbb{R}^n$; then,
the flow $\phi_{0,t}$ transports $D(0)$ to $D(t)=\phi_{t,0}(D(0))$.
A problem of interest is to understand how the volume
$v(t)=\operatorname{vol}(D(t))$ evolves with $\phi_{t,0}$.
$v(t)$ satisfies the equation
$$\begin{align}
\begin{matrix}
\dot{v}(t)&=&
\int\limits_{D(0)}(\nabla_y\cdot
f)(t,\phi(t;0,y))
\det\left[\frac{\partial\phi}{\partial y}(t;0,y)\right] dy\\
&=&\int\limits_{D(t)}\nabla_y\cdot f(t,y) dy
\end{matrix}
\tag{4}\label{liouvfor}
\end{align}
$$
To see this, apply the change of variables formula for integration,
to obtain
$$
\begin{align}
\begin{matrix}
v(t)&=&\int_{D(t)} dy=\int_{\phi_{t,0}D(0))} dy\\
&=& \int_{D(0)} \left|\det\left[\frac{\partial \phi_{t,0}}{\partial
y}(y)\right]\right| dy
\end{matrix}\tag{5}\label{multchange}
\end{align}
$$
Since $\phi_{t,0}(\cdot)$ is a family of diffeomorphisms and
$\phi_{0,0}=Id$, we can ignore the absolute value in \eqref{multchange}. Differentiating with respect to $t$ gives
$$
\dot{v}(t)= \int_{D(0)} \frac{d}{dt}
\det\left[\frac{\partial \phi_{t,0}}{\partial
y}(y)\right] dy
$$
The conclusion then follows from the statement of the OP along with another application of the change of variables formula.
