This paper presents an investigation of modeling and solving system of differential equations in the study of mechanical systems with holonomic constraints. A method is developed for constracting equation of motion for mechanical system with constraints. A technique is developed how to approximate the solution of the problem that is obtained from modeling of kinematic constraint equation which is stable. A perturbation analysis shows that velocity stabilization is the most efifcient projection with regard to improvement of the numerical integration. How frequently the numerical solution of the ordinary differential equation should be stabilized is discussed. A procedure is indicated to get approximate solution when the systems of differential equations can’t be solved analytically. A new approach is applied for constructing and stabilyzing Runge-Kutta numerical methods. The Runge-Kutta numerical methods are reformulated in a new approach. Not only the technique of formulation but also the test developed for its stability is new.Finally an example is presented not only to demonstrate how the stability of the solution depends on the variation of the factor but also how to find an approximate solution of the problem using numerical integration.