  We present four frequently used finite difference methods and establish the
error bounds for the discretization of the Dirac equation in the massless and
nonrelativistic regime, involving a small dimensionless parameter $0<
\varepsilon \ll 1$ inversely proportional to the speed of light. In the
massless and nonrelativistic regime, the solution exhibits rapid motion in
space and is highly oscillatory in time. Specifically, the wavelength of the
propagating waves in time is at $O(\varepsilon)$, while in space it is at
$O(1)$ with the wave speed at $O(\varepsilon^{-1}).$ We adopt one leap-frog,
two semi-implicit, and one conservative Crank-Nicolson finite difference
methods to numerically discretize the Dirac equation in one dimension and
establish rigorously the error estimates which depend explicitly on the time
step $\tau$, mesh size $h$, as well as the small parameter $\varepsilon$. The
error bounds indicate that, to obtain the `correct' numerical solution in the
massless and nonrelativistic regime, i.e. $0<\varepsilon\ll1$, all these finite
difference methods share the same $\varepsilon$-scalability as time step
$\tau=O(\varepsilon^{3/2})$ and mesh size $h=O(\varepsilon^{1/2})$. A large
number of numerical results are reported to verify the error estimates.
