This thesis presents a new approach for linearization of large multibody dynamic systems. The approach uses an analytical differentiation of terms evaluated in a numerical equation formulation. Because the method is based on a relative coordinate formalism, it is more efficient than any finite difference method without the concern of determining the proper dithering values. This new approach was generalized to include closed-loop systems, damping and eigenvalue sensitivities with respect to design variables. A number of examples were used to illustrate the accuracy and efficiency of the algorithm.