Polynomial with subset of critical points and values prescribed Motivated by this question I wish to pose the following question:
Given $k$ points $(x_1, y_1), \ldots (x_k, y_k)$ with (WLOG) $x_i < x_{i+1}$, can we find a polynomial $p(x)\in\Bbb R[x]$ satisfying


*

*$p(x_i) = y_i$

*$p'(x_i) = 0$

*For $x_1 < x < x_k$ we have $p'(x) = 0 \iff x\in \{x_1, \ldots, x_k\}$


Note the key difference with the linked problem in that we do not demand all of the critical points and values be prescribed, indeed the polynomial may have as many as we wish, and possibly a bunch of non-real roots. We simply demand that we know all of the real roots on some interval $[x_1, x_k]$. I am also willing to grant that $y_i=y_j\implies i = j$ carte blanche in the setup, as this was an issue in the original problem.
A "stretch goal" of sorts is that in an ideal world we would also have some control over the nature of the critical points, eg. prescribe a condition such as $p'(x_i)p'(x_{i+1})<0$ or possibly even greater control.

What does it all mean? This would answer one of the so-called "calculus superstition" problems undergraduates often have, i.e. any snippet of a random graph is secretly a polynomial graph, we teaching them are just making them suffer needlessly.

One approach is to start with $p_0(x)$ determined by the derivative criteria, using basic calculus to produce $P_0(x)$ satisfying our second condition, and choosing a primitive so that say $p(x_1)=y_1$. Now if we scale $p_0$ in addition, we can make it so $p(x_2) = y_2$ as well if needed, as we might as well consider $x_1=0$ by a variable shift, if necessary. But for $3$ or more points, this approach no longer works.
One is tempted to attempt the following construction in the case no $y_i =0$ (which I'm willing to grant for the time being as this won't even work in the special case):
Construct $P_0(x)$ as before, and consider $A(x)$ so that $A(x_i) = y_i / P_0(x_i)$, so that $P_1(x) = A(x_i)P_0(x_i) = y_i$. We immediately run into a problem though as here $P_1'(x) = A'(x)P_0(x) + A(x)P_0'(x)$ which makes this problem self-referential: we'd need to be able to construct $A(x)$ with the same critical points and it's own values there prescribed.
For my part I'm not sure I believe such a thing exists, but unlike the original problem the freedom in degree frees up a lot of wiggle room. The only answer to the linked question given (at the time of this posting) is Alexandre Eremenko's and of course those familiar with interpolation recognize its truth and that the two free parameters mean one could prescribe all of one and $2$ of the other with a given degree, but it's not clear to me (possibly because I am missing something quite obvious) that with no degree bound doesn't make the modified version of the problem work.
 A: I'm pretty sure that what you request does exist and sketch an approach. However I have some doubt the result would satisfy your goals. You want a polynomial which, on the interval, has exactly the critical points $(x_i,y_i)$ you desire and is strictly increasing or decreasing  between them. Presumably there is a continuous target function you wish to mimic. 
Note that $y_i=y_{i+1}$ would force a critical point in between but otherwise it is no problem to have $y_i=y_j.$ 
However one can imagine a polynomial which is increasing or decreasing as required, but nearly horizontal, in the long but thin rectangles with corners  $(x_{i-1}+\varepsilon, \frac{y_{i-1}+y_i}2\pm\varepsilon)$ and   $(x_{i}-\varepsilon, \frac{y_{i-1}+y_i}2\pm\varepsilon)$  and increasing/decreasing/flat as required, with the exact critical points, in the appropriate short tall rectangles with top (or bottom)  $(x_i,y_i \pm \varepsilon).$ The approach I sketch below might have similar behavior except with the abrupt nearly vertical transitions at the midpoints of sub-intervals. Perhaps also specifying the inflection points and signs of the second derivative (if possible) would be more satisfactory.
Let me first mention a few things which might work better. Hermite interpolation allows you to specify $f(x_i)$ and $f'(x_i)$ at $n$ points with a polynomial of degree $2n-1.$ This might well come out well, and one could add further points and/or higher order derivatives. However I don't know that there is a way to see this always works.
Using Bernstein Polynomials one can get a constructive proof that for any continuous function $f(x)$ on $[a,b]$ and $\varepsilon \gt 0$ there is a polynomial $p(x)$ (possibly of very high degree) with $|f(x)-p(x)| \lt \varepsilon$ on $[a,b].$ However this is not designed to hit specific points.
Here then is my approach (which is shakier than I thought but still seems sound): The desired polynomial can be $p(x)=y_0+\int_{x_0}^xq(t) dt$ where $q$ is a polynomial which is positive/negative/zero in the right places and such that $\int_{x_{i}}^{x_{i+1}}q(t)dt=y_{i+1}-y_{i}.$ The idea is to make $n-1$ polynomials where $q_i(x)$ is positive/negative/zero in the right places but nearly zero except on a very short sub-interval near the middle of $(x_i,x_{i+1})$  so the $(n-1) \times (n-1)$ matrix $M$ with  values $\int_{x_j}^{x_{j+1}}q_i(t)$ has $\pm 1$ on the diagonal and off diagonal entries very close to $0.$ Then $M^{-1}$ is a matrix with similar structure (diagonal entries near $\pm 1$ and small off diagonal entries) and can be used with the vector ${\mathbf y}$ of values $y_{i+1}-y_i$ to get the positive vector ${\mathbf c}=M^{-1}{\mathbf y}$ whose entries can be used for  $q(x)=\sum_1^{n-1}c_iq_i(x).$ 
You know that you want  $q$ (and the $q_i$) to be zero only in certain places. Using $\prod_1^n(x-x_i)$ achieves that. But you also want it to be positive (or negative) in between. This can be achieved by using $\pm\prod_1^n(x-x_i)^{e_i}$ where the sign in front is chosen appropriately and $e_i$ is $1$ where the sign changes and $2$ where it stays the same. Multiply this by some very small constant $c \gt 0$ to keep the $y$ values between $-0.1$ and $0.1$, perhaps raise it to a (large) odd power to further diminish  it. Call the result $q^*(x).$ Next consider the helper polynomial $h(x)=(1-x^2)^m$ where $m$ is a very large integer and the domain is $[-1,1].$ Then $h(\pm 1 )=0$, $h(1)=1$ and $h(x)$ is positive but very small outside of $(-\varepsilon,\varepsilon).$ Now let $q_i=q^{*}h_i$ where $h_i$ is a transformation of $h$ shifted and stretched so that the domain covers all of $[2x_1-x_n,2x_n-x_1]$ , the peak is at $x=(x_i+x_{i+1})/2$, and is just high enough to make $\int_{x_i}^{x_{i+1}}q_i(t)=\pm 1.$ If all the "very small" quantities are small enough this seems as if it would work, though not be the most satisfying or practical route.
