  We explore higher-dimensional generalizations of the
Runge-Kutta-Wentzel-Kramers-Brillouin method for integrating coupled systems of
first-order ordinary differential equations with highly oscillatory solutions.
Such methods could improve the performance and adaptability of the codes which
are used to compute numerical solutions to the Einstein-Boltzmann equations. We
test Magnus expansion-based methods on the Einstein-Boltzmann equations for a
simple universe model dominated by photons with a small amount of cold dark
matter. The Magnus expansion methods achieve an increase in run speed of about
50% compared to a standard Runge-Kutta integration method. A comparison of
approximate solutions derived from the Magnus expansion and the
Wentzel-Kramers-Brillouin (WKB) method implies the two are distinct
mathematical approaches. Simple Magnus expansion solutions show inferior long
range accuracy compared to WKB. However we also demonstrate how one can improve
on the standard Magnus approach to obtain a new "Jordan-Magnus" method. This
has a WKB-like performance on simple two-dimensional systems, although its
higher-dimensional generalization remains elusive.
