This paper presents a computer-based method for formulation and efficient solution of nonlinear, constrained differential equations of motion for spatial dynamic analysis of mechanical systems. Nonlinear holonomic constraint equations and differential equations of motion are written in terms of a maximal set of Cartesian generalized coordinates, three translational and four rotational coordinates for each rigid body in the system, where the rotational coordinates are the Euler parameters. Euler parameters, in contrast to Euler angles or any other set of three rotational generalized coordinates, have no critical singular cases. The maximal set of generalized coordinates facilitates the general formulation of constraints and forcing functions. A Gaussian elimination algorithm with full pivoting decomposes the constraint Jacobian matrix, identifies dependent variables, and constructs an influence coefficient matrix relating variations in dependent and indpendent variables. This information is employed to numerically construct a reduced system of differential equations of motion whose solution yields the total system dynamic response. A numerical integration algorithm with positive-error control, employing a predictor-corrector algorithm with variable order and step size, integrates for only the independent variables, yet effectively determines dependent variables.