A method is presented for formulating and solving the Newton-Euler equations of motion of a system of interconnected rigid bodies. The digital simulation may involve numerical integration of the kinematic equations as well as the dynamic equations. The reaction forces and torques resulting from rigid constraints imposed at the connecting joints are also determined. The derivation of kinematic expressions for first and higher derivatives is demonstrated based on direct differentiation of the rotation matrix in the spirit of the classical vector approach. A representative problem in spatial mechanism analysis is solved and illustrated with numerical results.