Using Generalized Inverse Method presented by Udwadia and Kalaba in 1992, we can obtain equations to exactly describe the motion of constrained systems. When the differential equations are numerically integrated by any numerical integration scheme, the numerical results are generally found to veer away from satisfying constraint equations. Thus, this paper deals with the numerical integration of the differential equations describing constrained systems. Based on Baumgarte method, we propose numerical methods for reducing the errors in the satisfaction of the constraints.