We discuss a novel strategy for computing the eigenvalues and eigenfunctions of the relativistic Dirac operator with a radially symmetric potential. The virtues of this strategy lie on the fact that it avoids completely the phenomenon of spectral pollution and it always provides two-side estimates for the eigenvalues with explicit error bounds on both eigenvalues and eigenfunctions. We also discuss convergence rates of the method as well as illustrate our results with various numerical experiments.