We are concerned with the numerical integration of ODEs of the form
$\epsilon^2 \psi_{xx} + a(x)\psi=0$ for given $a(x)\ge\alpha>0$ in the highly oscillatory regime $0<\epsilon\ll 1$ (appearing as a stationary Schr\"odinger equation, e.g.). In two steps we derive an accurate finite difference scheme that does not need to resolve each oscillation:
1) With a WKB-ansatz the dominant oscillations are "transformed out", yielding a much smoother ODE.
2) For the resulting oscillatory integrals we devise an asymptotic expansion both in $\eps$ and $h$. In contrast to existing strategies, the presented method has (even for a large spatial step size $h$) the same weak limit (in the classical limit $\epsilon\to 0$) as the continuous solution. Moreover, it has an error bound of the order $O(\epsilon^2 h^2)$.