  This paper is concerned with an efficient numerical method for solving the 1D
stationary Schr\"odinger equation in the highly oscillatory regime. Being a
hybrid, analytical-numerical approach it does not have to resolve each
oscillation, in contrast to standard schemes for ODEs. We build upon the
WKB-based marching method from [1] and extend it in two ways: By comparing the
$\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$ methods from [1] we design an
adaptive step size controller for the WKB method. While this WKB method is very
efficient in the highly oscillatory regime, it cannot be used close to turning
points. Hence, we introduce for such regions an automated method coupling,
choosing between the WKB method for the oscillatory region and a standard
Runge-Kutta-Fehlberg 4(5) method in smooth regions. A similar approach was
proposed recently in [2, 3], however, only for a $\mathcal{O}(h)$-method.
Hence, we compare our new strategy to their method on two examples (Airy
function on the spatial interval $[0,\,10^{8}]$ with one turning point at $x=0$
and a parabolic cylinder function having two turning points), and illustrate
the advantages of the new approach w.r.t.\ accuracy and efficiency.
