Abstract: In this talk we deal with the numerical solution of differential equations of the type ε²u″(x) + a(x)u(x) = 0 . Since for a small parameter ε > 0 the solution exhibits rapid oscillations with a wave length proportional to ε, standard schemes for ODEs become inefficient as they would have to resolve every oscillation. We present our main strategies for the construction of a third order (w.r.t. the step size) WKB-based numerical scheme: First, by using well known WKB-techniques from quantum mechanics, the given ODE can be reformulated as a less oscillatory one. Second, we need a special treatment for the highly oscillatory integrals, which appear in the Picard-approximation of the solution. We present global error bounds for the case of exactly as well as numerically computed phases. For an explicitly available phase our method has the property of asymptotical correctness w.r.t. the small parameter ε. We conclude with a discussion of some numerical examples.