In this article, we propose a new algorithm for solving the forward dynamics of multibody systems consisting of rigid bodies connected in arbitrary topologies by localized joints and/or soft links, possibly actuated or not. The simulation is based on the implicit time integration of the Lagrangian model of these systems, where the soft links are modeled by Cosserat rods parameterized by assumed strain modes. This choice imposes a predictor–corrector structure on the approach, and requires computing both the residual vector and the Jacobian of the residual vector of the dynamics constrained by the time integrator. These additional calculations are handled here with a new Newton–Euler recursive inverse dynamics algorithm and its linearized tangent version. The approach is illustrated with numerical examples from the Cosserat rod literature and from recent robotic applications.