Research into joint modelling methods has grown substantially over recent years. Previous research has predominantly concentrated on the joint modelling of a single longitudinal outcome and a single time-to-event outcome. In clinical practice, the data collected will often be more complex, featuring multiple longitudinal outcomes. Harnessing all available measurements in a single model is advantageous and should lead to improved and more specific model predictions. Notwithstanding the increased flexibility and better predictive capabilities, the extension of the classical univariate joint modelling framework to a multivariate setting introduces a number of technical and computational challenges. These include the high-dimensional numerical integration required, modelling of multivariate unbalanced data, and estimation of standard errors. Consequently, software capable of fitting joint models to multivariate data has been lacking. Building on recent methodological developments, we describe the extension of the classical joint model to multiple continuous longitudinal outcomes and describe how to fit it using a Monte Carlo Expectation-Maximization algorithm. Practical issues of Monte Carlo integration, convergence, and standard error estimation are discussed. The work is integrated into the R package, joineRML (https://CRAN.R-project.org/package=joineRML), which will also be described.