The Jacobian matrix of the equations of motion of a system using velocity transformation technique is derived via variation methods to apply the implicit integration algorithm, DASSL. The concept of generalized coordinate partitioning is used to parameterize the constraint set with independent generalized coordinates. DASSL is applied to determine independent generalized coordinates and velocities. Dependent generalized coordinates, velocities, accelerations and Lagrange multipliers are explicitly retained in the formulation to satisfy all of the governing kinematic and dynamic equations. The derived Jacobian matrix of a system is proved to be valid and accurate both analytically and through solution of numerical examples.