Find a smooth function with prescribed moments In several contexts I’ve encountered variants of the following problem :
let $m_0,m_1,m_2$ be real numbers such that $0 < m_1 < m_0$ and
$\frac{m_1^2}{m_0} <m_2 < m_1$. Then, show that there is a nonnegative smooth
function $f : [0,1] \to {\mathbb R}^{+}$ such that
$$
\int_{0}^{1} f(x)dx=m_0, \ 
\int_{0}^{1} xf(x)dx=m_1, \
\int_{0}^{1} x^2f(x)dx=m_2.
$$
The easiest version is when "smooth" means just "piecewise continuous" ; in that case
the only proof I know uses staged functions with complicated parameters. Is there
a proof that avoids such tedious computations ?
I guess that one can always find a ${\cal C}^{\infty}$ solution $f$, and even
an analytic one, but I have no clear idea on how to proceed. 
My (very vague) thoughts :
Bump functions ?
Integration by parts, transforming the problem into an easier interpolation problem ?
EDIT (09/24/2014) : In answer to Han de Bruijn’s comment, if one works on a general $[a,b]$ instead of $[0,1]$, so that $m_k=\int_a^b x^k f(x)dx$, then the inequalities
become
$$
am_0 < m_1 < bm_0, \ \frac{m_1^2}{m_0} < m_2 < (a+b)m_1-abm_0
$$
Note that the inequation $\frac{m_1^2}{m_0} < m_2$ (which remarkably
contains no $a$ or $b$) comes from the Cauchy-Schwarz inequality.
 A: One can find an analytic function [or even a polynomial] with prescribed moments. 
Let's normalize to $m_0=1$. Let $V$ be the set of nonnegative real-analytic functions [or just polynomials] $f:[0,1]\to\mathbb R^+$ with $\int_0^1 f=1$. The goal is to show that the image of $V$ under the map $$\Psi(f) = \left(\int_0^1 xf(x)\,dx, \int_0^1 x^2f(x)\,dx\right)$$
covers the region $D=\{(x,y):0<x<1,\ x^2<y<x\}$. 
Since $\Psi$ is a linear map, $\Psi(V)$ is a convex subset of $\mathbb R^2$. Let $\Gamma = \{(t,t^2):0\le t\le 1\}$ and observe that $D$ is the interior of the convex hull of $\Gamma$. By the Lemma stated below, it suffices to prove that $\Gamma\subset \overline{\Psi(V)}$.
To this end, take a sequence of nonnegative rational functions [or polynomials] that converge to Dirac delta $\delta_t$ in the sense of distributions. For example, 
$$
f_n(x) = \frac{c_n}{1+n(x-t)^2}
$$
where $c_n$ is chosen so that $\int_0^1 f_n=1$. [Replace $f_n$ with an approximating  polynomial, if desired.] One does not really need to talk about distributions here: just observe that for every $\delta>0$, 
$$
\int_{[0,1]\setminus [t-\delta,t+\delta]}f_n\to 0
$$
and conclude that 
$$
\int_0^1 x^k f_n(x)\,dx \to t^k,\quad k=1,2,\dots 
$$

Lemma
Let $A$ be a convex subset of $\mathbb R^n$. If $B\subset \overline{A}$, then the interior of the convex hull of $B$ is contained in $A$. 
Proof
Suppose $b\not\in A$ is an interior point of the convex hull of $B$. Since $A$ is convex,  there is a hyperplane $L$ through $b$ such that $A$ lies in a closed halfspace $H$ with $\partial H=L$. Since $B\subset \overline{A}\subset H$, the convex hull of $B$ is contained in $H$, hence cannot have  $b$ as an interior point. Contradiction.  
