In this paper, we study multiplicity of semi-classical solutions to nonlinear Dirac equations of space-dimension \begin{document}$ n $\end{document} : \begin{document}$ \begin{equation*} -i\hbar\sum\limits_{k = 1}^n \alpha_k \partial_k u+a\beta u+V(x)u = f(x,|u|)u,\; \text{in} \mathbb{R}^n, \end{equation*} $\end{document} where \begin{document}$ n\geq 2 $\end{document} , \begin{document}$ \hbar>0 $\end{document} is a small parameter, \begin{document}$ a>0 $\end{document} is a constant, and \begin{document}$ f $\end{document} describes the self-interaction which is either subcritical: \begin{document}$ W(x)|u|^{p-2} $\end{document} , or critical: \begin{document}$ W_{1}(x)|u|^{p-2}+W_{2}(x)|u|^{2^*-2} $\end{document} , with \begin{document}$ p\in (2,2^*), 2^* = \frac{2n}{n-1} $\end{document} . The number of solutions obtained depending on the ratio of \begin{document}$ \min V $\end{document} and \begin{document}$ \liminf\limits_{|x|\rightarrow \infty} V(x) $\end{document} , as well as \begin{document}$ \max W $\end{document} and \begin{document}$ \limsup\limits_{|x|\rightarrow \infty} W(x) $\end{document} for the subcritical case and \begin{document}$ \max W_{j} $\end{document} and \begin{document}$ \limsup\limits_{|x|\rightarrow \infty} W_{j}(x), j = 1,2, $\end{document} for the critical case.