  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.
