The validation and parallel implementation of a numerical method for the solution of the time-dependent Dirac equation is presented. This numerical method is based on a split operator scheme where the space-time dependence is computed in coordinate space using the method of characteristics. Thus, most of the steps in the splitting are calculated exactly, making for a very efficient and unconditionally stable method. We show that it is free from spurious solutions related to the fermion-doubling problem and that it can be parallelized very efficiently. We consider a few simple physical systems such as the time evolution of Gaussian wave packets and the Klein paradox.The numerical results obtained are compared to analytical formulas for the validation of the method.

3 In [13], µ 2 0 = 1.51 GeV 2 . This is the scale at which the gluon distribution is initialized, and it mostly effects the normalization of the energy density. 4 Qs/(g 2 µ) ≈ 0.75 in [13]