  A numerical method is developed to solve the time-dependent Dirac equation in
cylindrical coordinates for 3-D axisymmetric systems. The time evolution is
treated by a splitting scheme in coordinate space using alternate direction
iteration, while the wave function is discretized spatially on a uniform grid.
The longitudinal coordinate evolution is performed exactly by the method of
characteristics while the radial coordinates evolution uses Poisson's integral
solution, which allows to implement the radial symmetry of the wave function.
The latter is evaluated on a time staggered mesh by using Hermite polynomial
interpolation and by performing the integration analytically. The cylindrical
coordinate singularity problem at $r=0$ is circumvented by this method as the
integral is well-defined at the origin. The resulting scheme is reminiscent of
non-standard finite differences. In the last step of the splitting, the
remaining equation has a solution in terms of a time-ordered exponential, which
is approximated to a higher order than the time evolution scheme. We study the
time evolution of Gaussian wave packets, and we evaluate the eigenstates of
hydrogen-like systems by using a spectral method. We compare the numerical
results to analytical solutions to validate the method. In addition, we present
three-dimensional simulations of relativistic laser-matter interactions, using
the Dirac equation.
