We describe a method for calculating the roots of special functions
satisfying second order linear ordinary differential equations. It exploits the
recent observation that the solutions of a large class of such equations can be
represented via nonoscillatory phase functions, even in the high-frequency
regime. Our algorithm achieves near machine precision accuracy and the time
required to compute one root of a solution is independent of the frequency of
oscillations of that solution. Moreover, despite its great generality, our
approach is competitive with specialized, state-of-the-art methods for the
construction of Gaussian quadrature rules of large orders when it used in such
a capacity. The performance of the scheme is illustrated with several numerical
experiments and a Fortran implementation of our algorithm is available at the
author's website
