OK, here is my current thinking on Pietro's approach.

We use the differential equation
$$
(x^4-34x^3+x^2)U'''(x) + 3(2x^3-51x^2+x)U''(x)+(7x^2-112x+1)U'(x)+(x-5)U(x)=0
\tag{ODE3}$$
Multiply by $x^k$, $k \ge 0$, then integrate by parts as much as possible. The result is
$$
\int \!{x}^{k-1} \left( {k}^{3}{x}^{2}-34\,{k}^{3}x
+3\,{k}^{2}{x}^{2}+{k}^{3}-51\,{k}^{2}x+3\,k{x}^{2}-27\,kx+{x}^{2}-5\,
x \right) U \left( x \right) {dx}
=
\left( {x}^{k+2}{k}^{2}-34\,{x}^{1+k}{k}^{2}+{x}^{k+2}k+{x}^{k}{k}^{2
}-17\,{x}^{1+k}k+{x}^{k+2}-10\,{x}^{1+k} \right) U \left( x \right)
+
\left( -{x}^{1+k}k-{x}^{3+k}k+34\,{x}^{k+2}k+{x}^{1+k}+2\,{x}^{3+k}-
51\,{x}^{k+2} \right) U' \left( x \right)
+
\left( {x}^{k+4}-34\,{x}^{3+k}+{x}^{k+2} \right) U'' \left( x \right)
\tag{A}$$
In particular, for $k=0$,
$$
\int (x-5)U(x)\;dx
=
\left( {x}^{2}-10\,x \right) U \left( x \right) + \left( 2\,{x}^{3}-
51\,{x}^{2}+x \right) U' \left( x \right) + \left( {x}^{
4}-34\,{x}^{3}+{x}^{2} \right) U'' \left( x
\right)
\tag{A0}$$
Write $Q_k(x)$ for the right-hand-side of (A) and $\int R_k(x)U(x)\;dx$ for the left side. We will be doing differences like
$$
\big[Q_k(x)\big]_a^b := Q_k(b)-Q_k(a)
$$
since that will equal the integral on the left $\int_a^b R_k(x)U(x)\;dx$.

We want to arrange a solution $U$ so that $\int_0^c R_k(x)U(x)\;dx=0$ for $k \ge 0$. As Pietro noted, this will give us the recurrence we want for the moments $M(k):=\int_0^c x^kU(x)\;dx$. And in particular from (A0) we would have $\int_0^c (x-5)U(x)\;dx = 0$ so that $M(1)=5M(0)$.

Now consider (ODE3). First look at solutions as we approach $x=c$ from the left. Three linearly independent solutions are asymptotically: $1$, $(c-x)^{1/2}$, $(c-x)$. Other terms are higher powers (integer and half-integer). My calculations show:

for the term $1$, we get $Q_k(c^-)=1/2\, \left( 24+17\,\sqrt {2} \right) \left( 24\,k+24+7\,\sqrt {2} \right) \left( 17+12\,\sqrt {2} \right) ^{k}$

for the term $(c-x)^{1/2}$, we get $Q_k(c^-)=0$, amazingly

for the term $(c-x)$, we get $Q_k(c^-)=\left( -9792-6924\,\sqrt {2} \right) \left( 17+12\,\sqrt {2} \right) ^{k}$

terms $(c-x)^2$ and higher all produce $Q_k(c^-)=0$.

Now look at solutions as we approach $x=c_o$ from the right. Three linearly independent soltuions are asymptotically: $1$, $(x-c_o)^{1/2}$, $(x-c_o)$.

for the term $1$, we get $Q_k(c_o^+)=1/2\, \left( 17\,\sqrt {2}-24 \right) \left( -24\,k-24+7\,\sqrt {2} \right) \left( 17-12\,\sqrt {2} \right) ^{k}$

for the term $(x-c_o)^{1/2}$, we get $Q_k(c_o^+)=0$

for the term $(x-c_o)$, we get $Q_k(c_o^+)=-12\,\sqrt {2} \left( 17-12\,\sqrt {2} \right) ^{k+2}$

for higher terms, we get $Q_k(c_o^+)=0$

Now look at solutions as we approach $x=c_o$ from the left. Three linearly independent soltuions are asymptotically: $1$, $(c_o-x)^{1/2}$, $(x-c_o)$.

for the term $1$, we get $Q_k(c_o^-)=1/2\, \left( 17\,\sqrt {2}-24 \right) \left( -24\,k-24+7\,\sqrt {2} \right) \left( 17-12\,\sqrt {2} \right) ^{k}$

for the term $(c_o-x)^{1/2}$, we get $Q_k(c_o^-)=0$

for the term $(x-c_o)$, we get $Q_k(c_o^-)=-12\,\sqrt {2} \left( 17-12\,\sqrt {2} \right) ^{k+2}$

for higher terms, we get $Q_k(c_o^-)=0$

Now look at terms as we approach $x=0$ from the right. Three linearly independent terms are asymptotically: $1$, $\log x$, $(\log x)^2$.

for the term $1$, we get $Q_k(0^+)=0$

for the term $\log x$, we get $Q_k(0^+)=0$

for the term $(\log x)^2$, we get $Q_k(0^+)=0$ for $k \ge 1$ but $Q_0(0^+)=2$

for higher terms (such as $x$, $x\log x$, etc.), we get $Q_k(0^+)=0$

After all that, we come to the strategy for a solution. Start with the solution of (ODE3) approaching $c$ from the left that is asymptotically $(c-x)^{1/2}$. If we use only that, then $Q_k(c^-)=0$. (At the end we will multiply by a constant to get our final result.) Follow this solution to the left to the next singularity, that is approaching $c_o$ from the right. It will look like $r_1 + r_2(x-c_o)^{1/2}+r_3(x-c_o)+\dots$ for some constants $r_1, r_2, r_3$. Now consider solutions approaching $0$ from the right. Using only the two basis elements asymptotically $1$ and $\log x$, following them to the right until we approach $c_o$ from the left, we can choose a linear combination so that we get $r_1$ and $r_3$ from before: $r_1+r_2^*(c_o-x)^{1/2}+r_3(x-c)+\dots$. We want the same $r_1$ and $r_3$ as from the right, but $r_2^*$ is probably not the same as $r_2$. This way we get $Q_k(c_o^+)=Q_k(c_o^-)$ and $Q_k(0^+)=0$. So that all boundary terms cancel as desired.

Now we still have a constant factor to set. If $\int_0^c U(x)\;dx$ turns out to be nonzero, choose that factor so that $\int_0^c U(x)\;dx = 1$. Then (as noted at the beginning) we will have $A_n = \int_0^c x^n U(x)\;dx$ for all $n \ge 0$. Another thing we have to hope for is that $U(x)$ turns out to have constant sign: The way we constructed it, we have $U(x)>0$ near $c$, but we have to hope it turns out to be nonnegative everywhere.

I will see if I can get Maple to produce some pictures. I'm doing fine from $c$ going left to $c_o$, and from $c_o$ going both ways, but Maple seems reluctant to do (ODE3) from $0$ going up.

**added**

I am no expert on numerical ODEs. But here is what I get. Starting at the right, using the solution of (ODE3) with $U(x) = (c-x)^{1/2}+O((c-x)^{3/2})$ as $x \to c^-$, when I reach $c_o$ it looks like
$$
U(x) = 1831.7 - 5769.99(x-c_o)^{1/2}-13916.8 (x-c_o)+O((x-c_o)^{3/2})
\qquad\text{as } x \to c_o^+
$$
Then, starting with the soution satisfying
$$
U(x) = 1831.7 + 0.000103 (c_o-x)^{1/2}-13916.8 (x-c_o)+O((c_o-x)^{3/2})
\qquad\text{as } x \to c_o^-
$$
we end up with $0$ on the $(\log x)^2$ term:
$$
U(x) = O(\;|\log x|\;)
\qquad\text{as } x \to 0^+
$$
The graph seems to be positive everywhere! So we merely need to divide by the integral.

Here is a plot. (Not very informative.)

Here it is near the point $c_o$:

9more comments