A new representation of solutions to the equation −y′′+q(x)y=ω2y is obtained. For every x the solution is represented as a Neumann series of Bessel functions depending on the spectral parameter ω. Due to the fact that the representation is obtained using the corresponding transmutation operator, a partial sum of the series approximates the solution uniformly with respect to ω which makes it especially convenient for the approximate solution of spectral problems. The numerical method based on the proposed approach allows one to compute large sets of eigendata with a nondeteriorating accuracy.