  We consider continuous Dirac operators defined on $\mathbf{R}^d$,
$d\in\{1,2,3\}$, together with various discrete versions of them. Both
forward-backward and symmetric finite differences are used as approximations to
partial derivatives. We also allow a bounded, H\"older continuous, and
self-adjoint matrix-valued potential, which in the discrete setting is
evaluated on the mesh. Our main goal is to investigate whether the proposed
discrete models converge in norm resolvent sense to their continuous
counterparts, as the mesh size tends to zero and up to a natural embedding of
the discrete space into the continuous one. In dimension one we show that
forward-backward differences lead to norm resolvent convergence, while in
dimension two and three they do not. The same negative result holds in all
dimensions when symmetric differences are used. On the other hand, strong
resolvent convergence holds in all these cases. Nevertheless, and quite
remarkably, a rather simple but non-standard modification to the discrete
models, involving the mass term, ensures norm resolvent convergence in general.
