In this talk, I will review our recent works on numerical methods and analysis for solving the Dirac equation in the nonrelativistic limit regime,

involving a small dimensionless parameter which is inversely proportional to the  speed of light. In this regime, the solution is highly oscillating in time and the energy becomes unbounded and indefinite, which bring significant difficulty in analysis and heavy burden in numerical computation. We begin with four frequently used finite difference time domain (FDTD) methods and the time splitting Fourier pseudospectral (TSFP) method and obtain their rigorous error estimates  in the nonrelativistic limit regime by paying particularly attention to how error bounds depend explicitly on mesh size and time step as well as the small parameter.Then we consider a numerical method by using spectral method for spatial derivatives combined with an exponential wave integrator (EWI) in the Gautschi-type  for temporal derivatives to discretize the Dirac equation. Rigorious error estimates show that the EWI spectral method has much better temporal resolution than the FDTD methods for the Dirac equation in the nonrelativistic limit regime. These methods and results are then extended to the nonlinear Dirac eqaution in the nonrelativistic limit regime. Finally, we present a multiscale time integrator Fourier pseudospectral (MTI-FP) method for the Dirac equation, which is uniformly accurate for the dimensionless paramter. Numerical results demonstrate that our error estimates are sharp and optimal. This is a joint work with Yongyong Cai, Xiaowei Jia, Qinglin Tang and Jia Yin.