Interpolate 4 points by an increasing polynomial I need to create a polynomial function that passes through the points $(0,0)$, $(x_1,y_1)$, $(x_2,y_2)$, $(x_3,y_3)$ where
$$ x_1 < x_2 < x_3\quad \text{  and } \quad y_1 < y_2 < y_3 $$
with the only constraint that that function has to be increasing for all $x > 0$. I don't mind the degree of the polynomial function.
I am using now a spline approach to build such a function, but I am trying to substitute by a polynomial. The Lagrange interpolation does not ensure that the function is increasing. 
 A: The task permits full solution in case of points $$(0,0), (x_1,y_1), (x_2,y_2),\quad 0<x_1<x_2,\quad 0<y_1<y_2.$$
Let $$t=\dfrac{x}{x_2},\quad f=\dfrac{y}{y_2},$$ then points choosing to  $$(0,0), (t,f), (1,1),\quad t,f\in[0,1].$$
Let consider polynomial function $$y(x)=ax+bx^2+cx^3,$$ $$y'(x)=a+2bx+3bx^2,\quad y'(x)>0, \text{ if }x\in (0,1).$$
Known that minimal value of derivation achieves on segment bounds or in inflection point if it belongs to segment $(0,1)$. Coordinates of this point can be calculated from equation $y''(t_i) = 0, \quad 2b+3bt_i = 0$ or $t_i = -\dfrac b{3c}$, so
$$\begin{cases}
at + bt^2 + ct^3 = f\\
a + b + c = 1\\
a > 0\\
a + 2b + 3c > 0\\
a - \dfrac{b^2}{3c} > 0, \text{ if }-\dfrac {b}{3c}\in(0,1),
\end{cases}$$
$$\begin{cases}
a + b + c = 1\\
b(t-t^2) + c(t-t^3) = t-f\\
b + c < 1\\
b + 2c > -1\\
b + c + \dfrac{b^2}{3c} < 1, \text{ if } bc\in(-3c^2,0)
\end{cases}.$$
The last inequality can be written as
$$
\left[
\genfrac{.}{.}{0pt}{0}
{\begin{cases}
c < 0\\
b \in (0, -3c)\\
b^2+3bc+3c^2-3c > 0
\end{cases}}
{\begin{cases}
c > 0\\
b\in(-3c,0)\\
b^2+3bc+3c^2-3c < 0
\end{cases}}
\right.$$
Discriminant of polinomial $b^2+3bc+3c^2-3c$ (when $c$ regarded as a parameter) is
$$d = (3c)^2 - 4(3c^2-3c) = 12c-3c^2.$$
$d<0$ when $c\not\in[0,4],$ so:
$$\left[
\genfrac{.}{.}{0pt}{0}
{\begin{cases}
c < 0\\
b \in (0, -3c)\\
\end{cases}}
{\begin{cases}
c\in[0,4]\\
b\in[b_1,b_2]\\
\end{cases}}
\right.
$$
where 
$$b_1 = \max(-3c, -1.5c-\sqrt{3c-0.75c^2},$$
$$b_2 = \min(0, -1.5c+\sqrt{3c-0.75c^2}.$$
Built using Wolfram Alpha (1) graph shows that when $c\geq0$ this restriction It does not break the connection area of the allowed parameter values:

When $c<0,$ similar graph shows us allowed zone as two triangles with vertices in the points $(0,0), (0,1), (1.5,0.5)$ and $(0,-1),(0,0),(-0.5,0):$

General picture shows us allowed zone as two triangles with vertices in the points $(0,0), (0,1), (1.5,0.5)$ and $(0,-1),(0,0),(-0.5,0):$

If do not take into account the limitations of $c\in[0,4]$, the shape of the permitted area resembles forked spear.
Thus
$$\begin{cases}
a + b + c = 1\\
b = \dfrac{t-f}{t-t^2} - c(1+t)\\
b + c < 1\\
b + 2c > -1\\
\left[
\genfrac{.}{.}{0pt}{0}
{\begin{cases}
c < 0\\
b \in (0, -3c)\\
\end{cases}}
{\begin{cases}
c\in[0,4]\\
b\in[b_1,b_2]\\
\end{cases}}
\right.
\end{cases}.$$
Area restrictions defined, and it remains to determine whether to always find common points of the second equation of the system with this area. Since the second equation is infinite straight line, the only reason for the lack of a decision may be the output of the slope of this line beyond (-2,-1) angular coefficients of the second and third inequalities.
The form of the equation that this requirement is satisfied. Consequently, in the problem there are always values $b,c$ on the allowed area and $a=1-b-c$ in which the problem in this formulation has a solution.
A: Here are som thoughts abouth this question, but they do not lead to a conclusive statement. 
I will use the following evidence:

If a polynomial $P$ is such that $P(a)=b$, then $P$ is necessarily of the form
   $$ P(X)=(X-a)Q(X)+b,$$
  where $Q$ is a polynomial.

The idea is to apply this iteratively to $P$ such that $P(0)=0$, $P(x_i)=y_i$ with $1\leq i\leq3$. Since $P(0)=0$, $P$ is of the form
$$P=XA(X)$$
where $A$ is a polynomial.
Let us compute $A(x_1)$. We have $P(x_1)=x_1A(x_1)=y_1$, so we have 
$$ A(x_1)=\frac{y_1}{x_1}\equiv z_1.$$
We conclude that $A$ is of the form
$$A(X)=z_1+(X-x_1)B(X)$$
where $B$ is a polynomial. We can compute $B(x_2)$ using  $$y_2=P(x_2)=x_2A(x_2)=x_2\big[z_1+(x_2-x_1)B(x_2)\big].$$
We get $$B(x_2)=\frac{1}{x_2-x_1}\left(\frac{y_2}{x_2}-\frac{y_1}{x_1}\right)\equiv z_2.$$
We can likewise write $B$ as
$$B(X)=z_2+(X-x_2)C(X),$$
where $C$ is a polynomial. We compute $C(x_3)$ using
$$\begin{split}y_3&=P(x_3)=x_3A(x_3)=x_3\Big[z_1+(x_3-x_1)B(x_3)\Big]\\&=x_3\Bigg\{z_1+(x_3-x_1)\Big[z_2+(x_3-x_2)C(x_3)\Big]\Bigg\}.\end{split}$$
We obtain $$C(x_3)=\frac{1}{x_3-x_2}
\left[\frac{1}{x_3-x_1}\left(\frac{y_3}{x_3}-\frac{y_2}{x_1}\right)-z_2\right]\equiv z_3.$$
Finally (or at last) we can write $C$ as 
$$C(X)=z_3+(X-x_3)Q(X).$$ 
By necessary conditions, any polynomial satisfying the requirements
$P(0)=0$, $P(x_i)=y_i$ is of the form
$$P(X)=X\Bigg\{z_1+(X-x_1)\Big[z_2+(X-x_2)\big(z_3+(X-x_3)Q(X)\big)\Big]\Bigg\}.$$
Clearly if $Q=0$ then we have the Lagrange interpolation result.
If this polynomial does not work, then we must consider $Q\neq0$ with a positive leading coefficient. 
We must have $P'(0)\geq0$. Let us compute $P'(0)$. We have
$$P'(0)=A(0)=z_1-x_1z_2+x_1x_2z_3-x_1x_2x_3Q(0).$$
Consequently, we want
$$Q(0)\leq\frac{z_1-x_1z_2+x_1x_2z_3}{x_1x_2x_3}\equiv q_0.$$
If $q_0\leq0$ then necessarily $Q$ must be of degree at least 1.
By considering $P'(x_i)\geq0$, we can find an inequality for $Q(x_i)$.
This sums up to four inequalities, so the required minimal degree of $Q$
for the general case is probably at least $3$, unless I'm missing something.
I believe the general case requires a computer's help. 
This does not provide a full answer, I'm afraid. I anyway hope this could help. 
