The two-stage and three-stage Radau IIA stiff integrators, belonging to the implicit Runge–Kutta family, are implemented in a computationally efficient manner to solve flexible multibody systems. These problems feature large displacements and finite rotations together with small elastic deformations. The two-stage and three-stage algorithms are modified to integrate the finite element based, second order differential equations of multibody dynamics directly. A new simplified Newton iteration is implemented for the two-stage algorithm to reduce its computational cost. The resulting linear system of equations obtained at each Newton iteration is solved efficiently for both the two-stage and three-stage algorithms. The Jacobian matrix in the simplified Newton iterations is evaluated through numerical differentiation. A compact storage strategy is used to archive the banded, sparse matrices. The error estimation based on the embedded formula and step size selection are discussed in detail. The proposed schemes are validated with the help of different numerical examples. The simulation results are consistent with those of other types of integrators. These examples show that Radau IIA algorithms can solve stiff problems with acceptable speed while guaranteeing stability and accuracy.