An adaptive spectral method for oscillatory second-order linear ODEs with frequency-independent cost

We introduce an efficient numerical method for second order linear ODEs whose solution may vary between highly oscillatory and slowly changing over the solution interval. In oscillatory regions the solution is generated via a nonoscillatory phase function that obeys the nonlinear Riccati equation. We propose a defect-correction iteration that gives an asymptotic series for such a phase function; this is numerically approximated on a Chebyshev grid with a small number of nodes. For analytic coefficients we prove that each iteration, up to a certain maximum number, reduces the residual by a factor of order of the local frequency. The algorithm adapts both the step size and the choice of method, switching to a conventional spectral collocation method away from oscillatory regions. In numerical experiments we find that our proposal outperforms other state-of-the-art oscillatory solvers, most significantly at low-to-intermediate frequencies and at low tolerances, where it may use up to $10^6$ times fewer function evaluations. Even in high frequency regimes, our implementation is on average 10 times faster than other specialized solvers.