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.



Here are som thoughts abouth this question, but they do not lead to a conclusive statement.
I will use the following evidence:
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.