  A new representation for a regular solution of the perturbed Bessel equation
of the form $Lu=-u"+\left( \frac{l(l+1)}{x^2}+q(x)\right)u=\omega^2u$ is
obtained. The solution is represented as a Neumann series of Bessel functions
uniformly convergent with respect to $\omega$. For the coefficients of the
series explicit direct formulas are obtained in terms of the systems of
recursive integrals arising in the spectral parameter power series (SPPS)
method, as well as convenient for numerical computation recurrent integration
formulas.
  The result is based on application of several ideas from the classical
transmutation (transformation) operator theory, recently discovered mapping
properties of the transmutation operators involved and a Fourier-Legendre
series expansion of the transmutation kernel. For convergence rate estimates,
asymptotic formulas, a Paley-Wiener theorem and some results from constructive
approximation theory were used.
  We show that the analytical representation obtained among other possible
applications offers a simple and efficient numerical method able to compute
large sets of eigendata with a nondeteriorating accuracy.
