For a given integer $n\geq1$ and a bounded interval $I\subset(1,\infty)$ the goal is to determine all monic polynomials $P\in\mathbb{Z}[X]$ of degree $n$ with $n-1$ zeros in the interval $(0,1)$ and one zero in the interval $I$.
I have currently tried two main approaches: studying derivatives, and randomized sign checks.
Derivatives
The derivatives method works for $I=(1,N)$. Write $P(X)=X^n+c_{n-1}X^{n-1}+c_{n-2}X^{n-2}+\ldots+c_0$. Note that the zeros of the $(i+1)$'st derivative of $P$ lie between the zeros of the $i$'th derivative of $P$. Furthermore, note that the zero of $P^{(n-1)}(X)=n!X+(n-1)!c_{n-1}$ is $-c_{n-1}/n$, which is the arithmetic mean of the zeros of $P$, and thus at least the geometric mean $|c_0|^{1/n}\geq1$. It follows that any $i$'th derivative of $P$ has $i-1$ zeros in $(0,1)$ and one in $(1,N)$.
The algorithm is best explained through an example, so consider if $n=3$, so $P(X)=X^3-aX^2+bX-c$. Then $P''(X)=6X-2a$ needs a zero in $(1,N)$, giving $3<a<3N$. For any such $a$, $P'(X)=3X^2-2aX+b$ needs one zero in $(0,1)$ and one in $(1,N)$, giving $P'(0)>0$ and $P'(1)<0$, giving $0<b<2a-3$. Finally, for any such $b$, $P(X)$ needs two zeros in $(0,1)$ and one in $(1,N)$. If $\phi$ is the root of $P'(X)$ in $(0,1)$, we get $P(0)<0$, $P(\phi)>0$ and $P(1)<0$, giving two lower bounds and one upper bound on $c$.
So in general, for $i$ from $n-1$ to $0$, you iterate all possible values of $c_i$, and you keep track of all zeros of $P^{(i)}(X)$ in $(0,1)$ for the next value of $i$. With a depth first search, this finds all solutions. Note that this algorithm also proves there are always finitely many solutions.
This algorithm works great for low degrees, but quickly becomes terribly slow for high degrees, due to hitting many dead ends. I was able to find the polynomial $P(X)=X^6-66X^5+163X^4-155X^3+69X^2-14X+1$ with only real zeros, all in $(0,1)$ except for $63.47\ldots$ which is the smallest possible such value with degree $6$. For reference, the discriminant of $P$ is $3486377$. However I was unable to find a single example of degree $7$ with this algorithm.
Randomized
The idea behind this algorithm is very straight forward. Let $I=(a,b)$ for $1<a<b$. If $P$ suffices, then there exist values $0=x_0<x_1<\ldots<x_{n-1}=1$ on which $P$ alternates in sign, such that $P(1)$ and $P(a)$ are negative and $P(b)$ is positive. For fixed such values, looking for such polynomials $P$ is an integer linear programming problem. So we simply uniformly generate random such values and solve the ILP in hopes of finding solutions.
For $P$ with zeros $0<x_1<\ldots<x_{n-1}<1$ the probability of finding $P$ is equal to $$(n-2)!(x_2-x_1)(x_3-x_2)\ldots(x_{n-1}-x_{n-2})\geq(n-2)!\Delta_P^{1/2}(x_n-x_1)^{-1}(x_n-x_2)^{-1}\ldots(x_n-x_{n-1})^{-1},$$ which is at least $(n-2)!\Delta_P^{1/2}/x_n^{n-1}$. Using $\Delta_P\geq1$ thus gives us a very conservative bound.
By this algorithm, a potential candidate for smallest degree $7$ solution was found, namely $P(X)=X^7 - 154X^6 + 461X^5 - 554X^4 + 337X^3 - 108X^2 + 17X - 1$ with discriminant $394853489$ and all zeros in $(0,1)$ except for $150.97\ldots$, but due to the random nature, this is not necessarily the smallest possible. Note that the this polynomial has a relatively feasible $0.00288\ldots$ probability of being found, much better than the conservative estimate of $5!\cdot394853489^{1/2}/150.97^6=0.000000201\ldots$, in fact $5!(394853489^{1/2}/150.97^6)^{1/2}=0.00492\ldots$ is a much better estimate.
Newton's identities
One more idea considered but not implemented is to use Newton's identities. Let \begin{align} P(X)&=X^n-e_1X^{n-1}+e_2X^{n-2}-\ldots+(-1)^ne_0=(X-x_1)\ldots(X-x_n), \\e_k&=\sum_{i_1<\ldots<i_k}x_{i_1}\ldots x_{i_k}, \\p_k&=x_1^k+\ldots+x_n^k. \end{align} Then we have \begin{align} e_1&=p_1 \\2e_2&=e_1p_1-p_2 \\3e_e&=e_2p_1-e_1p_2+p_3 \\\vdots \end{align} So the coefficients of $P$ can be determined from the values $p_k$, and for them to be integers, each $p_k$ needs to satisfy an equation modulo $k$.
Then we know that $q_k=x_1^k+\ldots+x_{n-1}^k\in(0,n-1)$ for all $k$, and $p_k=q_k+x_n^k$, so $p_k\in(x_n^k,x_n^k+n-1)$ and $x_n^k\in(p_k-n+1,p_k)$. If $p_1,\ldots,p_{k-1}$ are chosen, then we know $p_k\in(x_n^k,x_n^k+n-1)\subset((p_{k-1}-n+1)^{k/(k-1)},p_{k-1}^{k/(k-1)}+n-1)$. If $I=(N,N+1)$, then $p_{k-1}$ is of the order $N^{k-1}$, which gives of the order $\frac{n}{k}N$ options for $p_k$. With $n$ options for $p_1$, this method thus tries of the order $\frac{n^n}{n!}N^{n-1}$ options. The most interesting consequence of this method is that this is thus also an upper bound on the total number of solutions.
Concrete question
This problem feels like there should be some literature on it, but I was unable to find any. More generally, I would expect there to be literature on determining all monic polynomials with its zeros is certain bounded intervals, specifically because this number is always finite. All the algorithms I used thus far are all right, but I would really like to look for higher degree polynomials, which is just infeasible with these methods.
As a final remark, consider writing $$P(X)=d_0X^n+d_1X^{n-1}(1-X)+\ldots+d_n(1-X)^n.$$ Then the coefficients appear to usually become much smaller. For example $P(X)=X^7 - 154X^6 + 461X^5 - 554X^4 + 337X^3 - 108X^2 + 17X - 1$ has $(d_0,\ldots,d_n)=(-1,9,-25,19,17,-27,10,-1)$. This seems to be related to Bernstein polynomials, due to all the zeros of $P$ in $(0,1)$.
If we consider $Q(X)=d_0+d_1X+\ldots+d_nX^n$ then $P(X)=Q(1/X-1)X^n$, so $Q$ needs to have $n-1$ zeros in $(0,\infty)$ and one zero in $(b^{-1}-1,a^{-1}-1)\subset(-1,0)$ for $I=(a,b)\subset(1,\infty)$. Furthermore, $P$ monic translates to $Q(-1)=1$. Maybe this translation somehow helps developing a quicker algorithm, due to the surprisingly small coefficients of $Q$.