  We apply the two-scale formulation approach to propose uniformly accurate
(UA) schemes for solving the nonlinear Dirac equation in the nonrelativistic
limit regime. The nonlinear Dirac equation involves two small scales
$\varepsilon$ and $\varepsilon^2$ with $\varepsilon\to0$ in the nonrelativistic
limit regime. The small parameter causes high oscillations in time which brings
severe numerical burden for classical numerical methods. We transform our
original problem as a two-scale formulation and present a general strategy to
tackle a class of highly oscillatory problems involving the two small scales
$\varepsilon$ and $\varepsilon^2$. Suitable initial data for the two-scale
formulation is derived to bound the time derivatives of the augmented solution.
Numerical schemes with uniform (with respect to $\varepsilon\in (0,1]$)
spectral accuracy in space and uniform first order or second order accuracy in
time are proposed. Numerical experiments are done to confirm the UA property.
