In this study, Constraint Force Equation Methodology (CFEM) is used to develop a multi-body dynamic analysis program with constraints. Seven constraint models are implemented to analyze constraint motions of multiple bodies. The augmented equations with the constraints are solved with the 4th order Runge-Kutta method for higher degree of accuracy. The analysis code is verified by comparing the analysis results of the motion of bodies with various constraints to published results.