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.