A fast and weakly stable method for computing the zeros of a particular class of hypergeometric polynomials is presented. The studied hypergeometric polynomials satisfy a higher order differential equation and generalize Laguerre polynomials. The theoretical study of the asymptotic distribution of the spectrum of these polynomials is an active research topic. In this article we do not contribute to the theory, but provide a practical method to contribute to further and better understanding of the asymptotic behavior. The polynomials under consideration fit into the class of Sobolev orthogonal polynomials, satisfying a four--term recurrence relation. This allows computing the roots via a generalized eigenvalue problem. After condition enhancing similarity transformations, the problem is transformed into the computation of the eigenvalues of a comrade matrix, which is a symmetric tridiagonal modified by a rank--one matrix. The eigenvalues are then retrieved by relying on an existing structured rank based fast algorithm. Numerical examples are reported studying the accuracy, stability and conforming the efficiency for various parameter settings of the proposed approach.
Mastronardi et al. (Wed,) studied this question.