  A special purpose solver, based on the Magnus expansion, well suited for the
integration of the linear three neutrino oscillations equations in matter is
proposed. The computations are speeded up to two orders of magnitude with
respect to a general numerical integrator, a fact that could smooth the way for
massive numerical integration concomitant with experimental data analyses.
Detailed illustrations about numerical procedure and computer time costs are
provided.
