We develop a fast, recursive and parameterization-free formulation for forward and inverse dynamics of soft robots modelled as multi-body systems with rigid and flexible bodies. The system is composed of bodies serially connected with discrete single-degree-of-freedom joints from a fixed base. The joint kinematics is captured through the exponential map of the Special Euclidean group SE(3). We couple the Newton-Euler equation for rigid bodies and a set of Partial Differential Equations (PDEs) on the SE(3) for dynamic Cosserat rods to formulate the dynamics of multi-body systems. Our proposed inverse dynamics recursively determines the system’s response as well as the required joint torques to follow a prescribed joint-space trajectory; while the forward dynamics inherits the same structure and determines the system’s motion for given joint torques. Unlike rigid multi-body systems, where the kinematic and dynamic recursions are decoupled, the inclusion of flexible bodies necessitates solving a coupled set of PDEs with applied separated Boundary Conditions (BCs). We develop a shooting-method-based BC solver to move BCs to one point and a numerical framework to integrate these equations. The framework adopts a fast finite difference method to semi-discretize the PDEs and transform them into ordinary differential equations, which are integrated using a geometrically exact integrator on the SE(3). Subsequent to the validation against multiple packages, we implement the proposed algorithms and study the dynamic response of an example. A software library, named SimUlator for Rigid-Flexible Robots (SuRFR), to simulate generic rigid-flexible multi-body systems, including the presented example, is distributed.