Abstract: With the development of science and technology, it is more and more important to study the dynamics of variable mass system such as jet aircraft and rocket, and it is always hoped that the solutions of the variable mass system are stable or asymptotically stable. It is difficult to study the stability by using Lyapunov direct methods because of the difficulty of constructing Lyapunov functions directly from the differential equations of the mechanical system. This paper presents an indirect method for studying stability, that is, gradient system method. This method can not only reveal the internal structure of dynamic system, but also help to explore the dynamic behavior such as the stability, asymptotic and bifurcation. The function V of the gradient system is usually taken as a Lyapunov function, so the gradient system is more suitable to be studied with the Lyapunov function. The equations of motion for the holonomic mechanical system with variable mass are listed, and all generalized accelerations are obtained in the case of non-singular system. A class of gradient system with negative-definite matrix is proposed, and the stability of the solutions of the gradient system is studied. This kind of gradient system and variable mass mechanical system are combined, then the conditions under which the solutions of the mechanical systems with variable mass can be stable or asymptotically stable are given. Further the mechanical system with variable mass whose solution is stable or asymptotically stable is constructed by using the gradient system with non-symmetrical negative-definite matrix. Through specific examples, it is studied that the solutions of the single degree of freedom motion of a variable mass system are stable or asymptotically stable under some conditions of the laws of mass change, particle separation velocity and force. The method is also suitable for the study of other constrained mechanical systems.