A method to calculate static solutions for mechanical systems containing rigid and flexible bodies modeled by finite elements is described. The formulation of the equations makes use of generalized strains, which leads to an extended set of equations for both these generalized strains and nodal coordinates, together with constraint equations imposing the relations between these two groups of coordinates. The associated Lagrangian multipliers are the generalized stresses. The resulting iteration scheme appears to be quite robust in comparison with more traditional methods, especially if some displacements are prescribed. Once a static solution has been found, the linearized equations of motion about this solution can be obtained in terms of a set of minimal coordinates, that is, in the degrees-of-freedom (DOFs). In addition, a continuation method is described for tracing a branch of static solutions if some parameters are varied. The method is of the familiar predictor–corrector type with a linear or cubic predictor and a corrector with a step size constraint. Applications to a large-deflection problem of a curved cantilever beam, large deflections of a fluid-conveying tube and its resulting instability, and the buckling of an overconstrained parallel leaf-spring mechanism due to misalignment are given.