  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.
