Is there some description of polynomial differential operators, $\mathcal{D}=\sum f_i(x) D_x^i$ such that, if $h$ is a polynomial all of whose roots are in $[0,1]$, then so are all the roots of $\mathcal{D} h$?
Some examples:
$\bullet$ $\mathcal{D} h = f(x) h(x)$, where $f$ is a polynomial all of whose roots are in $[0,1]$.
$\bullet$ $\mathcal{D} = D_x$, and more generally $D_x^k$, by Rolle's theorem.
$\bullet$ $\mathcal{D} h = (x-r)^{1-a} D_x \left( (x-r)^a h \right)$ for $r \in [0,1]$ and $0 \leq a$. By Rolle's theorem applied to $(x-r)^a h$. Properly, this only makes sense when $a$ is of the form $m/(2n+1)$, so that we can raise negative numbers to the $a$, but formally it gives a differential operator for all $a$ and it preserves real-rootedness by continuity.
$\bullet$ $\mathcal{D} h = (x-r)^{k-a} D_x^k \left( (x-r)^a h \right)$ for $r \in [0,1]$ and $k-1 < a$, generalizing the above.
$\bullet$ $\mathcal{D} h = \left( \prod_{i=1}^s (x-r_i)^{k-a_i} \right) \cdot D_x^k \left( \prod_{i=1}^s (x-r_i)^{\vphantom{k}a_i} h \right)$ for $0 \leq r_1 \leq r_2 \leq \cdots \leq r_s \leq 1$ and $a_1$, $a_2$, ..., $a_s \geq k-1$, generalizing the above, if I didn't make any errors.
$\bullet$ Any composition of the above.
Motivation: I just wrote an answer which came down to manipulating a differential operator until I could show that it had this property. (The question was about polynomials in $x^{-1}$ all of whose roots are in $[2, \infty)$, but that's easily equivalent.) I'm curious whether there is an algorithm for this sort of thing.
