Computing the Zeros and Turning Points of Solutions of Second Order Homogeneous Linear ODEs

Algorithms to compute the zeros and turning points of the solutions of second order ODEs $y^{\prime\prime}+B(x)y^{\prime}+A (x)y=0$ are developed. Two fixed point methods are introduced. The first method consists of fixed point iterations that are built from the first order differential system relating the problem function with a contrast function w; the contrast function w has zeros interlaced with those of y, and it is a solution of a different second order ODE. This method, which generalizes previous findings [J. Segura, SIAM J. Numer. Anal., 40 (2002), pp. 114--133], requires the evaluation of the ratio of functions y/w. The second method is based on fixed point iterations stemming from the second order ODE; it requires the computation of the logarithmic derivative $y^{\prime}/y$. Both are quadratically convergent methods; error bounds are provided. The particular case of second order ODEs depending on one parameter, $y_n^{\prime\prime}+B_n (x) y_n^{\prime}+A_n (x) y_n (x)=0$, with applications to the computation of the zeros and turning points of special functions, is discussed in detail. The combination of both methods provides algorithms for the efficient computation of the zeros and turning points of a broad family of special functions, including hypergeometric and confluent hypergeometric functions of real parameters and variables (Jacobi, Laguerre, and Hermite polynomials are particular cases), Bessel, Airy, Coulomb, and conical functions, among others. We provide numerical examples showing the efficiency of the methods.