Question: Given a polynomial ODE $\dot{x}=f(x)\in\mathbb{R}^n$ that possesses a Darboux polynomial$^*$ $p(x)$ satisfying $\dot{p}(x)=c(x)p(x)$ for some function $c(x)$ (called the cofactor) how can one write $f(x)$ in terms of $p(x)$ and/or $c(x)$?
Example: To make things a little more concrete with an example, let $H(x)$ be a first integral satisfying $\dot{H}(x)=0$ (this is the $c(x)=0$ case). Then an ODE possesses a first integral $H(x)$ iff it can be written in the following skew-gradient form $$\dot{x}=S(x)\nabla H(x),$$ for some matrix $S(x)=-S(x)^T$. My question is: is there an equivalent statement for the case of ODE possessing a Darboux polynomial?
My initial thoughts: you can let $$f(x) = \frac{c(x)p(x)}{\nabla p(x)^T A \nabla p(x)} A \nabla p(x)$$ for some matrix $A$, then $f(x)$ has $p(x)$ as a Darboux polynomial. However, this relies on $ \frac{c(x)p(x)}{\nabla p(x)^T A \nabla p(x)}A$ being polynomial for $f(x)$ to also be polynomial. Also this doesn't really fit with the above first-integral case when $c(x)\rightarrow 0$.
$^*$ Darboux polynoimals are also known as second integrals or weak integral for non-polynomial ODEs.
Edit: As suggested, including a skew symmetric matrix $S=-S^T$ would make $f$ satisfy the first integral condition on its level set ($c(x)\rightarrow 0$ or $p(x)\rightarrow 0$), that is:
$$f(x) = \left(\frac{c(x)p(x)}{\nabla p(x)^T A \nabla p(x)} A + S\right)\nabla p(x)$$
The following should solve your problem under the assumption that $\nabla p(x)$ is an everywhere non-vanishing vector field:
You want $\nabla p\cdot f = c p$, so let $\alpha = c p / |\nabla p|^2$, then $\nabla p \cdot (f - \alpha \nabla p) = 0$ and you are back to your first example. Find a skew symmetric $S$ so that $X:=f-\alpha \nabla p = S \nabla p$. Then $f=(\alpha I + S) \nabla p$ has the required form.
If the data in your question are given, then $f$ is already assumed to be a polynomial.
In order to obtain a skew symmetric map $S$ as above, at a point $x$ we define the linear map $S_x:{\Bbb R}^n \to {\Bbb R}^n$ by declaring (omitting the base point $x$ in the notation):
$$ S Y = \frac{X (\nabla p \cdot Y)- \nabla p (X \cdot Y)}{|\nabla p|^2} \ \ Y \in {\Bbb R}^n .$$ By construction, $Y \cdot (S_x Y)=0$ for all $x$ and every $Y\in {\Bbb R}^n$, so $S_x$ is represented by a skew-symmetric matrix in any orthonormal basis of ${\Bbb R}^n$. Also, as required, $S_x(\nabla p) = X$ since $X\cdot \nabla p = 0$.