The objective of this study is to present an accurate and simple method to describe the motion of constrained mechanical or structural systems. The proposed method is an elimination method to require less effort in computing Moore-Penrose inverse matrix than the generalized inverse method provided by Udwadia and Kalaba. Considering that the results by numerical integration of the derived second-order differential equation to describe constrained motion veer away the constrained trajectories, this study presents a numerical integration scheme to obtain more accurate results. Applications of holonomically or nonholonomically constrained systems illustrate the validity and effectiveness of the proposed method.