When large relative displacements take place between the bodies in a multi-body Wave Energy Conversion system linearisation of the constraints on motion imposed by the joints between the bodies is no longer valid and a non-linear timedomain analysis is necessary. As a part of the Techno-Economic Optimisation of Wave Energy Conversion (TEOWEC) software, which has been developed at the Centre for Ocean Energy Research (COER), NUI Maynooth, we developed an algorithm for the dynamic simulation of Multi-Body Systems for Wave Energy Conversion (MBS4WEC) with fully non-linear representation of the body-to-body joints. The algorithm is based on the Jointcoordinate formulation, which provides a systematic procedure to transform the mixed differential-algebraic equations of motion in body coordinates, for open chain systems, to a minimal set of ODEs. When a closed-loop chain occurs, the same method can be adopted by removing one or more kinematic joints from each loop. Knowing the topology of the system, a path matrix is generated and together with the formulation of data structures representing the body-to-body joints, the Velocity Transformation Matrix is computed. The main advantage of this approach is a fast and automatic generation of the Velocity Transformation Matrix, which leads to a higher computational efficiency, especially for complex systems. This paper presents the equations underpinning the method together with results for simulation of two specimen floating multi-body systems. These two are a simple multi-body hinged barge and a device with a sliding internal reaction mass. In each case the results are contrasted to the results produced by a linearised analysis of the same system.