  We introduce a double power series method for finding approximate analytical
solutions for systems of differential equations commonly found in cosmological
perturbation theory. The method was set out, in a non-cosmological context, by
Feshchenko, Shkil' and Nikolenko (FSN) in 1966, and is applicable to cases
where perturbations are on sub-horizon scales. The FSN method is essentially an
extension of the well known Wentzel-Kramers-Brillouin (WKB) method for finding
approximate analytical solutions for ordinary differential equations. The FSN
method we use is applicable well beyond perturbation theory to solve systems of
ordinary differential equations, linear in the derivatives, that also depend on
a small parameter, which here we take to be related to the inverse wave-number.
  We use the FSN method to find new approximate oscillating solutions in linear
order cosmological perturbation theory for a flat radiation-matter universe.
Together with this model's well known growing and decaying M\'esz\'aros
solutions, these oscillating modes provide a complete set of sub-horizon
approximations for the metric potential, radiation and matter perturbations.
Comparison with numerical solutions of the perturbation equations shows that
our approximations can be made accurate to within a typical error of 1%, or
better. We also set out a heuristic method for error estimation. A Mathematica
notebook which implements the double power series method is made available
online.
