I have asked several questions on math.SE in order to compute numerically the poles of high-degree Padé approximations for $e^{-x}$, because a computation directly from the polynomial coefficients has poor numerical stability. I finally blundered upon an approach alluded to a paper by Campos and Calderón and a related item mentioned in a paper by Grünbaum which states that for polynomials with simple roots
If $$P(x)=(x-x_1)(x-x_2)\ldots (x-x_n)$$ then at any root $x_k$ of $P$ we have $$P''(x_k) = 2P'(x_k)\sum\limits_{j \neq k}\frac{1}{x_k-x_j}$$
In other words,
$$\frac{P''(x_k)}{2P'(x_k)} = \sum\limits_{j \neq k}\frac{1}{x_k-x_j}$$
(I called this a "Stieltjes residue sum" in the title since it seems like his name pops up in this arena, and Grünbaum's article cites some work he did in electrostatics applications with these sorts of sums. If there's a more accurate name, please correct me + I'll change the title.)
For certain polynomials $P(x)$ satisfying a nonlinear second-order differential equation, the left side becomes very simple. For example, for the Bessel polynomials $y_n(x)$ and $\theta_n(x)$ with roots $y_k$ and $z_k$, respectively,
$$\frac{1}{y_j}+\frac{1}{y_j^2} = \sum\limits_{k\neq j} \frac{1}{y_k - y_j}$$
$$-1-\frac{n}{z_j} = \sum\limits_{k\neq j} \frac{1}{z_k - z_j}$$
and for the denominator $D_{pq}(x) = {}_1F_1(-q;-p-q;x)$ of the Padé approximation of $e^{-x}$ with roots $x_k$,
$$\frac{x_k+p+q}{2x_k} = \sum\limits_{k\neq j}\frac{1}{x_k-x_j}$$
These can be turned into nonlinear vector equations $\mathbf F(\mathbf x) = \mathbf 0$ where $\mathbf x$ is a vector of roots, and solved iteratively by the Newton-Raphson method given an initial guess that is vaguely close to the exact values. I've done this for the Bessel polynomial $\theta_n(x)$ (described in detail on math.stackexchange.com) and for $D_{pq}(x)$ and got very good results with fast convergence (8-10 iterations) even for 500th-degree polynomials.
What I don't understand is why this works so well. Can anyone explain why this form of a polynomial would converge so well? Finding roots of the polynomials given the coefficients is a very ill-conditioned problem. And this method has $n!$ possible solutions, so it seems like any iterative approach would be fraught with instability.
