The following properties define the polynomial $p(x)$ uniquely:
$\text{deg}(p(x))=7\\p(-1)=y_1,\ p'(-1)=d_{1,1},\ p''(-1)=d_{2,1},\ p'''(-1)=d_{3,1}\\ p(1)=y_2,\ \ \ \ p'(1)=d_{1,2},\ \ \ \ p''(1)=d_{2,2},\ \ \ \ p'''(1)=d_{3,2}$
The parameters $y_1,\ d_{1,1},\ d_{2,1},\ d_{3,1},\ y_2,\ d_{1,2},\ d_{2,2},\ d_{3,2}$ are real and chosen in such a way that the implication $-1\leq x\leq 1\Rightarrow-1\leq p(x)\leq1$ is true.
For each parameter, I need to find all real numbers that, when replacing the parameter, make the implication $-1\leq x\leq 1\Rightarrow-1\leq p(x)\leq1$ true.
There is always one solution, because, for instance, if I replace $y_1$ with a real number $r$, one solution is $r=y_1$. Moreover, the complete solution is an interval, $a\leq r\leq b$, because all points with $x\in[-1,1]$ will move the same way, when a parameter is modified monotonously.
How should this be done numerically? (I don't think there is an general exact solution)
I am not sure I understand the problem. The polynomial contains 8 coefficients and you give 8 conditions; so we have 8 linear equations for 8 unknowns. The system seems to have a unique solution; using a CAS, I have been able to generate the expressions for each coefficient as a function of the imposed conditions. These expressions are not complex at all; they are just linear combinations of the values of the imposed conditions.
Be sure I apologize not to put them here; I am almost blind and typing the formulas would be very hard for me. Please, tell me if I can do something else for helping you.
Happy New Year