WKB-based scheme with adaptive step size control for the Schr\"odinger equation in the highly oscillatory regime

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 (named after the physicists Wentzel, Kramers, Brillouin) marching method from [2] and extend it in two ways: By comparing the $\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$ methods from [2] 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 methods switching, 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 [9, 4], however, only for an $\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 on a parabolic cylinder function having two turning points), and illustrate the advantages of the new approach w.r.t.\ accuracy and efficiency.


Introduction
This paper deals with the numerical solution of the highly oscillatory 1D Schrödinger equation (1.1) Here, 0 < ε ≪ 1 is the rescaled Planck constant (ε := √ 2m ) and ϕ the (possibly complex valued) Schrödinger wave function. The real valued coefficient function a(x) is related to the potential. We shall assume here that it is bounded away from zero, i.e. a(x) ≥ τ for some τ > 0. The (local) wave length of a solution ϕ to (1.1) is given by λ(x) = 2πε √ a(x) . Hence, for small values of ε the solution is highly oscillatory, especially in the semi-classical limit ε → 0. Oscillatory problems like (1.1) appear in a wide range of applications, e.g., quantum mechanics, electron transport in semiconductor devices, and acoustic scattering. For instance, the state of an electron that is injected with the prescribed energy E from the right boundary into an electronic device (e.g., diode), modeled on the interval [0, 1], is described by the equation (see [2]) where k(x) := ε −1 E − V (x) is the wave vector and V denotes the electrostatic potential. Note that our assumption a(x) ≥ τ > 0 implies E > V (x), so the solution ψ E becomes oscillatory. Then, one is often interested in macroscopic quantities like the electron density n and the current density j, which are given by where f represents the injection statistics of the electrons. Here, ℑ(·) denotes the imaginary part and E(k) = ε 2 k 2 + V means the energy for a given wave vector k. In order to compute the quantities (1.3), the Schrödinger equation (1.2) has to be solved many times, as a fine grid in E(k) is needed. Hence, efficient methods for the solution of (1.2) are of great interest in such applications. Instead of solving the boundary value problem (1.2) directly, one can also solve the equivalent initial value problem, which results if equation (1.1) on the interval (0, 1) is augmented with the initial values ϕ(0) = ϕ 0 = 1 and εϕ ′ (0) = ϕ ′ 0 = − i a(0) with a(x) = E − V (x). The solution ϕ of this initial value problem and that one of problem (1.2) are then related by according to [2]. Indeed, the method proposed in this paper will deal only with initial value problems for the Schrödinger equation (1.1), but through this equivalence it is equally suitable for solving problem (1.2).
In [2,13,14], efficient and accurate WKB-based (named after the physicists Wentzel, Kramers, Brillouin; cf. [15]) numerical schemes have been developed for (1.1) in the oscillatory regime. By transforming out the dominant oscillations, they allow to compute a solution using a coarse spatial grid with step size h > λ. In fact, the grid limitation can there be reduced to at least h = O( √ ε). Now, this work adds on top of the algorithm from [2] an adaptive step size control as well as a switching mechanism. This allows the algorithm to switch to a standard ODE method (e.g., Runge-Kutta) during the computation in order to avoid technical or efficiency problems in regions where the coefficient function a(x) is very small or indeed equal to zero. We recall that the WKB-approximation is not valid close to turning points, i.e. where a(x) = 0. A switching mechanism was also used in [9], where the authors presented another WKB-based numerical scheme for the initial value problem corresponding to (1.1). Therefore, one goal of this work is to compare numerical results from our method with the results given by the method from [9], by considering two examples where the analytical solution is known. Since numerical methods for oscillatory problems is an active field of research, let us mention some references that are intended for more general oscillatory problems, and hence also include adequate methods for (1.1): first the two monographs [18,10]. The adiabatic integrators of [10,§XIV] are in fact closely related to a zeroth order WKB-approximation, see (2.2)-(2.3), below. Concerning (highly oscillatory) Hamiltonian boundary value methods we cite [7,1]. This paper is organized as follows: In Section 2 we give a short review of the second order (w.r.t. ε) WKB-marching method from [2]. Section 3 then describes the adaptive step size control algorithm as well as the switching mechanism. In Section 4 we recap the Runge-Kutta-WKB method from [9] and point out the difference between their step size control and the one used in this paper. In Section 5 we present numerical investigations on the error estimators of our algorithm as well as a comparison of the numerical results of our method and the method from [9]. We then conclude in Section 6.

The WKB-marching method
We aim at solving the Schrödinger equation (1.1), augmented with the initial conditions with some x 0 ∈ R. First we shall review the basics of the second order (w.r.t. ε) WKB-marching method from [2] with focus on the algorithm. The motivation for this method was the construction of a numerical scheme that is uniformly correct in ε and sometimes even asymptotically correct, i.e. the numerical error goes to zero with ε → 0 while the grid size h remains constant. For further details see [2]. The method consists of two parts: 1. Analytic pre-processing of (1.1) by transforming the equation into a smoother (i.e. less oscillatory) problem that can be solved accurately and efficiently on a coarse grid. 2. Obtaining a numerical scheme by discretization of the smoother problem.
Analytic pre-processing. The well-known WKB-approximation (cf. [15]) for the oscillatory regime where a(x) ≥ τ for some τ > 0, is based on inserting the ansatz into equation (1.1). After a comparison of ε-powers one obtains the first three functions φ p (x) as Here the symbols ± and ∓ in (2.3) and (2.5) correspond to the fact that there is always a pair of approximate solutions to the Schrödinger equation (1.1), by analogy to the two fundamental solutions of (1.1). Hence the general solution is then a linear combination of the two. Therefore, a truncation of the sum (2.2) after p = 2 leads to the second order (w.r.t. ε) asymptotic WKB-approximation of ϕ(x): with constants c 1 , c 2 ∈ C to be determined from initial or boundary conditions, and the phase is In the WKB-marching method of [2] this second order WKB-approximation is used to transform (1.1) into a smoother problem. To clarify our terminology, we point out that this method (as well as the Runge-Kutta-WKB method of Section 4) has both a WKB-order (w.r.t. ε; referring to the used cut-off in the asymptotic expansion (2.2)) and a numerical order (w.r.t. the step size h; referring to the convergence order). Firstly, using the notation the second order differential equation (1.1) with the initial conditions (2.1) can be reformulated as a system of first order differential equations: Here, the two matrices A 0 and A 1 are given by Then, the first order system (2.9) for U(x) is transformed by the change of variables with the two matrices where φ ε is the phase function defined in (2.7). This leads to the system 11) where N ε (x) is a (Hermitian) matrix with only off-diagonal non-zero entries: Since the transformation (2.10) eliminated the dominant oscillations, the system (2.11) can be solved numerically on a coarse grid {x n , n ∈ N 0 }. Then the original solution can be recovered by the inverse transformation It should be noted that, throughout the whole transformation from U(x) to Z(x), the phase integral φ ε (x) is assumed to be known exactly. For a generalization using a spectral method to numerically compute the phase see [5].
Numerical scheme. The derivation of the second order (in h) scheme for (2.11) is obtained via the second order Picard approximation (2.12) Since the entries of N ε (x) are highly oscillatory, (2.12) involves (iterated) oscillatory integrals. With φ ε assumed to be known exactly, they are then approximated using similar techniques as the asymptotic method in [12]. The first order (in h) scheme for (2.11) is derived by only taking into account the first two terms from (2.12). For both of these schemes we introduce the following notations: Further, let {x 0 , x 1 , ..., x N } be a grid we want to compute the solution on, and h := max 1≤n≤N |x n − x n−1 | be the step size. Then both schemes read as follows: First order scheme: Let Z 0 := Z I be the initial condition and let n = 0, ..., N− 1. Then the algorithm updates as with the (Hermitian) matrix and the phase increments Second order scheme: Let Z 0 := Z I be the initial condition and let n = 0, ..., N − 1. Then the algorithm updates as with the (Hermitian) matrix and the diagonal matrix

Step size control and switching mechanism
The WKB scheme is efficient in the highly oscillatory regime, but not applicable close to turning points, i.e. zeros of a(x), see [3]. This is evident already from the transformation (2.8), which does not make sense when a(x) ≤ 0. For mixed problems, e.g., the Airy equation on R + 0 (see Section 5.1), which has a turning point at x = 0, it is therefore convenient to couple two different methods: a method for highly oscillatory ODEs (e.g., WKBbased) away from the turning point, and a standard ODE method (e.g., Runge-Kutta) close to the turning point, where the solution is smooth anyhow. Here, we choose the well-known Runge-Kutta-Fehlberg 4(5) (RKF45) scheme (cf. [11]) as the standard ODE method. The latter method will be applied directly on equation (1.1) and not on (2.11), since the WKBtransformation (2.8)-(2.11) is not permitted at turning points. The exact switching mechanism as well as the introduction of an adaptive step size control to the algorithm will be described in the two following subsections.

The adaptive step size controller
In order to compute the solution efficiently, an adaptive step size controller, based on an estimator for the local truncation error, will be added to the numerical methods. This control allows the step size to increase or decrease while aiming to keep the error estimator as close as possible within a given error tolerance. To illustrate the functionality of this step size controller, we shall consider a numerical scheme of order k. We are then able to apply this step size control individually to the different methods mentioned above.
) be two numerical solutions of order k and k + 1 to approximate the exact solution Y (x n ) = (ϕ(x n ), ϕ ′ (x n )) of the initial value problem (1.1), (2.1). E.g., one could choose the WKB schemes (2.13) and (2.14) of h-order 1 and 2, respectively. Next we want to decide whether to accept the numerical solution at x n (typically the more accurate solution Y (k+1) n ) or rather to retry the computation with a modified step size. To this end we define the estimator for the local truncation error as Let h n,trial := x n − x n−1 be the (trial) step size which was used to compute the solutions at the current step n. We then use the common approach of varying the step size via the multiplicative control Here, we choose the factor θ n to be based on the so-called elementary controller (e.g., see [17]). Additionally, we introduce limitations in such way that the step size controller responds "smoothly" to abrupt changes in the solution behaviour, that is, the ratio between two consecutive step sizes should not be exorbitantly large or small. That said, we choose the factor similar to [6, p. 310] as where ATol = η·Tol and RTol = Tol are absolute and relative error tolerances for a given master tolerance Tol, the values 0.5 and 2 are design parameters that limit the ratio of two consecutive step sizes from below and above, and 0.9 is a common safety factor for increasing the probability of the next step to be accepted. Here, η is a scaling factor representing the gradual switchover between absolute and relative errors, depending on the behaviour of the solution. This is because for Y (k+1) n ∞ → 0 the ATol term in the numerator in (3.2) is dominating, whereas for large values of Y (k+1) n ∞ the RTol term is dominating. For our numerical simulations in Section 5 we choose η = 10 −2 . If ATol + RTol · Y (k+1) n ∞ ≥ est n , we accept the n-th step with the step size defined as h n := h n,trial and define the trial step size for the next step as However, if ATol + RTol · Y (k+1) n ∞ < est n , the n-th step gets rejected and a reattempt is done with the smaller (trial) step size h n,trial by updating its value as Since such an acceptance criterion is based on aiming to maintain the local error in each step as close as possible to the given error tolerances it is often referred to as error per step (EPS) control. In practice, using EPS control one can hope to achieve a global truncation error proportional to Tol k/(k+1) (e.g., see [6, p. 311]).

The switching mechanism
As already mentioned above, the algorithm shall automatically switch between two numerical methods, the WKB method in the oscillatory regime and another method, which is valid close to turning points. To realize this dynamical switching mechanism we follow a similar strategy as in [4]. To illustrate this procedure we now consider two numerical schemes of order k (1) and k (2) , where the superscripts (1) and (2) correspond to the two schemes. In each step, the adaptive step size algorithm from the previous section is applied to both schemes individually up to the definition of the quantities θ (1) n and θ (2) n , i.e., we just evaluate (3.1) and (3.2). Then, after checking the acceptance criterion ATol + RTol · Y n for each scheme (i) ∈ {(1), (2)}, the switching mechanism intervenes and it selects the acceptable numerical method that yields the larger value of θ (i) n , for i = 1, 2, hence yielding the larger proposed new step size. We thus favor the method with the smaller error estimator, discounted by its respective order k (i) . More precisely, we define and store the information on the method of choice in that n-th step. Through this procedure the algorithm does not only use the error estimator to find the next step size, but also to decide between the two methods. If at least one method was accepted, the algorithm sets Otherwise a reattempt is done with the smaller (trial) step size h n,trial by updating its value as h n,trial → Θ n · h n,trial .
We remark that the coupling of two methods, as presented above, could incur additional computational costs, since in every step both methods have to be applied in order to compute Θ n . However, in our case the WKB method (for highly oscillatory regions) and the standard ODE solver (for smoother regions) complement each other very well, yielding better results concerning accuracy as well as overall efficiency (see also 11 and 16).

The Runge-Kutta-WKB method
In this section we give a short review of the Runge-Kutta-WKB (RKWKB) method presented in [9] for the initial value problem (1.1), (2.1). The method is based on a dynamic switching mechanism between a standard Runge-Kutta scheme and a new stepping procedure that uses the WKB-ansatz (2.2) as an approximation of the true solution. This stepping procedure reads as follows: where and ϕ ′′ n is computed from ϕ n and equation (1.1). Here, f ± are chosen as WKB-approximations (2.2) of some finite order. For instance, one gets the second (WKB-)order method by setting Note that this choice is equal to (2.6), but one can easily choose f ± with higher (WKB-)orders. However, the stepping procedure (4.1)-(4.5) is always a first order (in h) numerical method. This holds because the coefficients γ ± and δ ± are chosen in such way that one finds from (4.1)-(4.3), as stated in [9]. It is also worth noting that in [9] the authors use a slightly different dynamic switching mechanism and a different step size control compared to the algorithm presented in Section 3. Firstly, their estimator of the relative error within the WKB-stepping procedure uses the difference between two numerical solutions ϕ (k) n and ϕ (k+1) n of different WKB-orders instead of different h-orders, simply since they do not have schemes of two different h-orders at their disposal. The algorithm then decides between a WKB step and a RK step by choosing the method with the smaller error estimator. In their more recent paper [4] the authors use a similar step size control and switching mechanism as presented in Section 3, but they do not limit the ratio of two consecutive step sizes and provide no option for controlling both absolute and relative errors as done in (3.2).
The goal of the following section is to compare numerical results from the WKB-marching method to results one gets using the RKWKB method. To both methods we will apply (exactly) the step size control and switching mechanism from Section 3, for the sake of comparability. Since the WKBstepping procedure of the RKWKB method is always of first order w.r.t. h, the definition of the error estimator (3.1) does not make sense. Therefore we shall use two different WKB orders (instead of h-orders) to be able to compute the error estimator in this case, as also done in [9]. Further, since our algorithm consists of a different step size update formula and switching criterion, we will call this modified RKWKB method simply RKWKBmod. Since this modification may produce slightly different numerical results compared to [9,4] we will also include the original RKWKB method from [9] into our comparison.

Numerical results: WKB-marching method vs. the Runge-Kutta-WKB method
In this section we will compare numerical results of the WKB-marching method and of the RKWKB method by applying both algorithms to two examples, where exact analytical solutions are available. The first example corresponds to a linear coefficient function a(x) and is taken from [9,4], whereas the second example involves a quadratic function a(x) and appears in [3]. In both examples the phase integral (2.7) in the WKB basis functions (2.3)-(2.5) can easily be computed exactly, hence we do not need any numerical integration routine here. By contrast, in [4] they evaluate the WKB basis functions with a numerical integration routine and therefore get another source of error in their method. Moreover, we will always use the second order WKB-marching method, since no scheme of higher WKB-order has been developed yet, see [2].
For clarity of the presentation we shall use in the sequel the following terminology for the methods to be compared: WKB+RKF45 WKB-marching method (see Section 2) + step size and switching algorithms to RKF45 (see Section 3) RKWKB original method from [9]: Runge-Kutta-WKB (see Section 4) + original step size and switching algorithms (see Section 4) RKWKBmod modified method from [9]: Runge-Kutta-WKB (see Section 4) + modified step size and switching algorithms (see Section 3) RKF45 Runge-Kutta-Fehlberg 4(5) scheme + step size algorithm (see Section 3) Table 1: Terminology for the methods to be compared.

First example: Airy function
The first example investigated in [9] is the Airy equation which results if one chooses the coefficient function a(x) = x. Here, Γ denotes the gamma function. The exact solution to the problem (5.1) is given by where Ai and Bi denote the Airy functions of first and second kind, respectively. This example demonstrates very well the advantages of a WKB method, since the solution becomes more and more oscillatory for large values of x. While standard adaptive ODE methods, e.g., Runge-Kutta methods, would need to decrease the step size more and more, a WKB-based method does not have to resolve the individual oscillations. Actually, it even allows to increase the step size for large x and is therefore highly efficient. This can be seen, e.g., in Figures 4-6.
Before starting to discuss the numerical solution of (5.1) we first remark that the evaluation of an oscillatory function (even of trigonometric functions such as sin or cos) is numerically ill posed for very large values of x. This is a generic problem and not related to the numerical solution of (5.1), or to the choice of the numerical method: Remark 5.1. We consider to evaluate ϕ exact (x) for an argument x that is only known with finite accuracy, specified by machine eps. Then, using the lowest order expansions for ϕ exact and ϕ ′ exact from Appendix A, we have x 3/2 ε as x → ∞ , and the relative error of ϕ exact (x(1 + eps)) is asymptotically eps x 3/2 /ε. Considering double precision, e.g., with Matlab's eps ≈ 2.2 · 10 −16 and x = 50, this generic evaluation error is 7.8 · 10 −10 for ε = 10 −4 and 7.8 · 10 −11 for ε = 10 −3 . These unavoidable errors limit the achievable accuracy, and it will play a role in Figures 10-11 below.
As a first step we shall now test the reliability of our choice of error estimator (3.1) -but only for the WKB steps of WKB+RKF45 and RKWKBmod. Since we know the exact solution to this problem, we can compute the local truncation error in each step and are able to compare it to the error estimator. Moreover, we apply the adaptive step size control from Section 3 with the error tolerance Tol = 10 −5 and set ε = 1.  According to the results of Figures 1 and 2, the error estimator is in excellent agreement with the local truncation error (for this example). Hence, it seems to be an adequate choice. We also find a very good agreement by plotting the respective relative errors for one single step as functions of the step size h for a fixed starting point x 0 , see Figure 3. Again, the relative error is for both methods much smaller than one. Note that in the case of WKB+RKF45 the relative error goes to zero, if h tends to zero. Therefore, the estimator seems to be asymptotically correct.  Table 2. The asymptotic expressions are based on well-known expansions, which can be found in Appendix A. The order for truncating the expansions was set to K = 3.  Now, we will give numerical results for solving the initial value problem (5.1) for the Airy equation with WKB+RKF45 and RKWKBmod. We recall that the algorithm automatically chooses between RKF45 and the respective WKB steps. To illustrate this difference in the Figures 4-9 as well as Figures 13-14, we will mark RKF45 steps with red dots, WKB steps using WKB+RKF45 with blue squares, and WKB steps using the RKWKBmod method with green triangles. To exclude the turning point x = 0, we solve the Airy equation (5.1) in Figures 4 and 5 on [0.1, 50].  According to Figures 4 and 5, WKB+RKF45 seems to perform slightly better than the RKWKBmod method in this example, in matters of global error. But at the same time WKB+RKF45 needs significantly fewer steps than RKWKBmod. Also, within the algorithm using WKB+RKF45, the switch from RKF45 steps to WKB steps happens earlier as can be seen in Figure 4. In Figure 6, the relative global errors of both algorithms are compared on [0.1, 10 8 ]. We find that the algorithm using WKB+RKF45 made fewer steps (58 vs. 91) while producing a slightly lower global error at the same time. The almost identical grid spacing of both methods from x ≈ 300 onwards is due to limiting the quotient of two consecutive step sizes, imposed in both methods. Furthermore, Figures 7 and 8 show the ratio between the WKB-and RKerror estimators as well as the ratio of the proposed step sizes (by the WKB and RKF schemes) for each step of the above simulations.  For both WKB+RKF45 and RKWKBmod, these plots demonstrate well the superiority of the RK scheme in the less oscillatory regime, i.e. for x small, as the ratio of the error estimators is very large there. The switching points to the WKB schemes are well defined (again for both WKB methods) -due to the monotonous behaviour of both the ratio of error estimators and the ratio of proposed step sizes. Note also that the step size ratios are bounded from above and below because of the limitations in (3.2).
We also want to give a perception of how WKB+RKF45 can perform, if the phase integral (2.7) can not be evaluated exactly. For this purpose we use the Clenshaw-Curtis quadrature (cf. [8]) to approximate the phase (2.7). That spectral method was already used in combination with the WKBmarching method, see [5, §5]. In [5] they approximate the phase at first on Chebyshev collocation points throughout the whole interval and use barycentric interpolation for the ODE-grid points x n . This was possible since they worked only on the "small" interval [0, 1]. In contrast, we will instead approximate the phase (2.7) in each interval [x n , x n+1 ] individually, since we are dealing here with the "long" interval [0.1, 10 8 ]. Figure 9 gives a comparison of the global error for the Airy equation (5.1), when using the exact phase vs. a numerically computed phase. According to these results the relative ϕerrors from computing the phase (2.7) numerically become visible only from x ≈ 10 7 onwards. We will next compare the numerical results of four methods, namely, WKB+RKF45, RKWKBmod, RKWKB, and RKF45 on the Airy equation for several values of the parameter ε. For RKF45, we do not use any smaller value of ε than 10 −1 , as the CPU time gets exorbitantly large. In Figure 10 we compare the accuracy of each of the methods depending on Tol. We find that WKB+RKF45 outperforms the other methods, in matters of global errors, particularly for small values of ε and Tol. By Remark 5.1 the accuracy limit at x = 50 is 7.8 · 10 −10 for ε = 10 −4 , which seems to explain the lower bound of the errors in Figure 10. Note also that for RKF45 the error increases like O(ε −1 ) for smaller ε-values, while it decreases for the WKB methods.
In Figure 11 we give a work-precision diagram; for a fair comparison between the methods, points showing the same error should be compared. There is a big difference between WKB+RKF45 and RKWKB(mod), regarding the CPU time: For ε = 10 0 , 10 −1 , 10 −2 this difference is particularly significant for small errors (stemming from small prescribed tolerances). For ε = 10 −3 WKB+RKF45 already beats RKWKB(mod) for all data points. For ε = 10 −4 the error intervals of WKB+RKF45 and RKWKB(mod) do not overlap. But for the same CPU time WKB+RKF45 yields much more accurate results than RKWKB or RKWKBmod. Overall we conclude from Figure 11 that WKB+RKF45 outperforms RKWKBmod and RKWKB significantly. Using RKF45, the CPU time increases drastically for smaller ε-values. Figure 12 displays the respective number of steps needed in each method, as a function of the prescribed tolerance. Note that for WKB+RKF45 and RKWKBmod the number of steps is bounded from below. This is because of limiting the quotient of two consecutive step sizes and it is clearly visible for ε = 10 −3 , 10 −4 .

Second Example: Parabolic cylinder function
The second example is taken from [3] and includes a quadratic coefficient function a(x): .
As before, let us compare numerical results only for WKB+RKF45 and RKWKBmod at first. There are two turning points, namely at x = 0 and x = 2. Therefore, we expect the two methods to make RKF45 steps near the turning points and WKB steps between them. Numerical results for the specific choice ε = 2 −6 are presented in Figures 13-14. According to Figures 13 and 14 no significant difference between WKB+RKF45 and RKWKBmod can be observed for ε = 2 −6 , in matters of global error. But again, WKB+RKF45 needs significantly fewer steps than RKWKBmod. Within the algorithm using WKB+RKF45, the switch-over between RKF45 steps and WKB steps happens closer to the turning points.
We shall now present numerical results for the four methods WKB+RKF45, RKWKBmod, RKWKB, and RKF45 for several values of ε. In Figure 15, we compare the global errors of each method depending on Tol. Here, we observe only a small difference between WKB+RKF45 and RKWKBmod. Using RKWKB, less smooth error curves can bee seen, with large peaks for lower errors (i.e. lower tolerances). In contrast to every other method, WKB+RKF45 seems to produce quite ε-independent global error curves. For RKF45 one sees that the error again increases like O(ε −1 ) for smaller values of ε.
In the work-precision diagrams in Figure 16, we observe for WKB+RKF45 that the CPU times are quite independent of ε, whereas for RKWKB(mod) they grow with decreasing ε, particularly for small errors (stemming from small prescribed tolerances). Overall we conclude from Figure 16 that WKB+RKF45 outperforms RKWKBmod and RKWKB (particularly for small ε and small tolerances) while showing an ε-uniform behaviour at the same time. Note also that, for RKF45, the CPU time increases drastically for smaller values of ε. Figure 17 shows the respective number of steps needed in each method, as a function of the prescribed tolerance. Remark 5.3. The computation of the reference solution involves the evaluation of the PCF, which is not readily available in Matlab. But the PCF can be related to the Kummer confluent hypergeometric function 1 F 1 (see [16, §13]), which is available in Matlab as hypergeom(). However, for small parameters ε this evaluation is very time consuming. For Figures 15-17 we thus computed the reference solutions by solving the initial value problem (5.2) with Matlab's routine ode45() (for all ε) using a very small error tolerance.

Conclusion
We have introduced in this paper an extension to the WKB-marching method from [2] by including into the algorithm an adaptive step size controller as well as a switching mechanism. In numerical simulations based on two examples this method yielded smaller global errors (particularly for small tolerances and small ε-values) in comparison to the Runge-Kutta-WKB method from [9,4], an alternative WKB-based scheme. Our tests revealed that the efficiency gain is mostly due to the different WKB method used here, while the different step size controls (here vs. [9,4]) do not play a big role. Our switching mechanism ensures well defined switching points between WKB steps and Runge-Kutta-Fehlberg 4(5) steps for oscillatory and, respectively, smoother regions of the ODE-solution. Especially for the Airy equation on the large spatial interval [0.1, 10 8 ] the efficiency of the method is demonstrated very well, as the scheme skips millions of oscillations within one step, while staying accurate at the same time. There is also a Matlab program available in a GitHub repository 1 , which also offers the possibility to compute the phase (2.7) numerically, as done for Figure 9.