Let $x\_1 < x\_2 < \ldots < x\_n$ and $y\_1 < y\_2 < \ldots < y\_n$ be two sequences
of $n$ real numbers. It is well known that there are polynomials that "interpolate"
in that $f(x\_i)=y\_i$ for all $i$, and the Lagrange interpolating polynomial
even warrants a solution of degree $ < n$. Now, what happens if we want the
polynomial $f$ to be nondecreasing on the interval $[x\_0,x\_n]$ ? Is there always
a solution, and is there a bound on the degree also ? 



Finite interpolation by a nondecreasing polynomial
This problem has appeared before in literature and is now well understood, I guess. The general version is when you have no restriction on the $y\_i$'s and you ask for an interpolating polynomial that is monotone on each sub-interval $[x\_ix\_{i+1}]$. The first paper proving the existence of such a polynomial is: 


W.Wolibner, "Sur un polynom d'interpolation", Colloq. Math (2) 1951, 136-137


but it is a non-constructive proof, as it uses the Weierstrass approximation theorem much like the answer given by Harald Hanche-Olsen above. Another proof for the case $0=y\_0\le \cdots \le y\_n=1$ is given in "Polynomial Approximations to Finitely Oscillating Functions" by W.J. Kammerer (Theorem 4.1) and the non-constructive aspect of his proof is the use of uniform convergence of appropriate Bernstein polynomials. In "Piecewise monotone polynomial interpolation", S.W. Young proves the same theorem and makes the final remark that the existence of such monotone interpolating polynomial is in fact *equivalent* to the Weierstrass theorem. On the other hand Rubinstein has some papers devoted to proving the existence of interpolating polynomials which are increasing in all of $\mathbb R$.


The first paper which gives bounds on the degrees is, I think, 


E. Passow, L. Raymon, "The degree of piecewise monotone interpolation", which is [here](http://www.jstor.org/stable/2040274)


and an improvement is made in "Exact estimates for monotone interpolation" by G.L. Iliev.
Note that the bounds are in terms of 
$$A=\max \Delta y\_i=\max (y\_{i}-y\_{i-1}) \qquad B=\min \Delta y\_i \qquad C=\min \Delta x \_i$$


And no uniform bound exists.


