I want to find a set of functions (if they exist) $ g_{l}(x)$ that are recursively defined as $$g_{l-1}=-(l+1)\rho'g_{l}$$ and satisfy the following integral equation up to index $l=k$: $$ 0 = \int_{-\infty}^{+\infty} (x^k)^{(l)}{d \over dx}( g_{l+1}\rho' )\rho^{l+1}dx $$ $\rho=\rho(x)$ is a known (but not necessarily named) probability density function that decays rapidly as $|x|\rightarrow\infty$, and $ (x^k)^{(l)} $ is shorthand for the $l-$th derivative of $x^k$. An obvious approach would be to find $g_{l+1}$, such that $ {d \over dx}( g_{l+1}\rho' ) = 0$ : $$g_{l+1}={C_{l+1}\over\rho'}$$ This solution, however fails to satisfy the recursive definition.
For $l=k$ the second equation simplifies to: $$0 = \int_{-\infty}^{+\infty} {d \over dx}( g_{k+1}\rho' )\rho ^{k+1}dx $$ which can be integrated by parts: $$0 = g_{k+1}\rho'\rho^{k+1}|_{-\infty}^{+\infty} - (k+1)\int_{-\infty}^{+\infty} g_{k+1}\rho'\rho^k dx$$ The first term goes to zero because of $\rho$, therefore we end up with: $$ 0 = \int_{-\infty}^{+\infty} g_{k+1}{d \over dx}(\rho^{k+1}) dx $$ At this point i'm stuck, my original hope was to find a function at the end of this iteration that could be propagated back easily, but this seems not to be the case.
As a bottom line, to give you some context: my goal is to find a function $f = g_0$ such that the evolution of $\rho(t,x)$ in the pde $${\partial \rho \over \partial t} = -{\partial \over \partial x}(\rho f) $$ is moment-preserving up to the k-th moment: $$ \mu_k = \int_{-\infty}^{+\infty}x^k \rho dx $$
For $k = 1 $, suitable functions are $f(x) = C\rho' $. The sign of C determines the behaviour (C<0 leads to a broadening of the distribution over time (see figure 1), C>0 leads to a concentration of the distribution).

Thanks in advance!