Abstract: Dynamic equations of flexible multibody systems are usually a set of stiff differential equations. At present, the common numerical method for solving the stiff differential equations filters out the high frequency by using the numerical damping. The computational efficiency of this method is still unsatisfactory. In order to reduce the stiffness of dynamic equations of flexible multibody systems so greatly that the equations can be solved by regular ordinary differential equation (ODE) solvers such as MATLAB ODE45 solver, methods of filtering high frequency vibrations during the process of modeling are studied. Stresses of flexible bodies are homogenized by their mean value over a time interval from now to a short time later. The homogenized stress is then employed to replace its origin when computing the virtual deformation power. In this way, the obtained model of the flexible multibody system will not contain harmful high frequency elastic vibrations. The range of frequencies can be controlled by the length of the time interval used to homogenize stresses. As validated by the numerical examples in this paper, the precision and efficiency of the proposed method are comparable to some stiff ODE solvers. Moreover, it works well when the stiff ODE solver fails to give correct solutions in a reasonable time. Comparisons of numerical examples show that the proposed method can be a new available approach to numerical analysis of flexible multibody systems.