A representation of solutions of the one-dimensional Dirac equation is obtained. The solutions are represented as Neumann series of Bessel functions. The representations are shown to be uniformly convergent with respect to the spectral parameter. Explicit formulas for the coefficients are obtained via a system of recursive integrals. The result is based on the Fourier-Legendre series expansion of the transmutation kernel. An efficient numerical method for solving initial-value and spectral problems based on this approach is presented with a numerical example. The method can compute large sets of eigendata with non-deteriorating accuracy.

22 pages, 1 figure. Added Corollary 4.9 and Remarks 5.1, 5.2