Computational run times in distinct element method (DEM) simulations of granular flow can be large and limit the size of the system being modeled. Hence, it is important to use efficient numerical integration schemes. This paper investigates some numerical integration schemes for their accuracy, stability and computational efficiency. It also investigates the effect of different particle contact damping algorithms on the model mathematical accuracy and stability. It is shown that the half-step leapfrog Verlet algorithm is the best integration scheme, while Enter is poor in terms of accuracy. Non-linear damping has been shown in the literature to be more realistic in terms of experimental data on particle impact coefficient of restitution. This was reproduced here. This paper also shows that non-linear damping reduces the mathematical error in the integration scheme because the force change is less discontinuous. However, in particle assembly simulations, filling a hopper, the non-linear damping model was less stable, probably because less energy is dissipated at low velocities.