I want to simulate the behaviour of a 2-DOF robotic manipulator, which is described by the following model:

$M(q)\ddot{q}=-C(q,\dot{q})\dot{q}-G(q)+\tau \text{}\text{}\text{}\text{}\text{}(1)$

Considering the fact that the 2x2 mass matrix M is positive definite, I could use the inverse matrix and break down the problem into the 4 first order ordinary differential equations and simulate it:

${x}_{1}={q}_{1}\Rightarrow \dot{{x}_{1}}={x}_{2}$

${x}_{2}=\dot{{q}_{1}}\Rightarrow \dot{{x}_{2}}=-{M}^{-1}(1,:)\cdot C\cdot {\left[\begin{array}{cc}{x}_{2}& {x}_{4}\end{array}\right]}^{T}-{M}^{-1}(1,:)\cdot G+{M}^{-1}(1,:)\cdot \tau $${x}_{3}={q}_{2}\Rightarrow \dot{{x}_{3}}={x}_{4}$

${x}_{4}=\dot{{q}_{2}}\Rightarrow \dot{{x}_{4}}=-{M}^{-1}(2,:)\cdot C\cdot {\left[\begin{array}{cc}{x}_{2}& {x}_{4}\end{array}\right]}^{T}-{M}^{-1}(2,:)\cdot G+{M}^{-1}(2,:)\cdot \tau $

Suppose I would like to use a solver that takes as an argument the mass matrix (a MATLAB ODE solver in particular) and don't use its inverse because this will also simplify the computation of the jacobian (I intend to simulate a 7-DOF robotic manipulator after that so providing the mass matrix would be great). How can I write the initial equation (1) as a series of first order ordinary differential equations and be able to simulate its response by using some software solvers ?

$M(q)\ddot{q}=-C(q,\dot{q})\dot{q}-G(q)+\tau \text{}\text{}\text{}\text{}\text{}(1)$

Considering the fact that the 2x2 mass matrix M is positive definite, I could use the inverse matrix and break down the problem into the 4 first order ordinary differential equations and simulate it:

${x}_{1}={q}_{1}\Rightarrow \dot{{x}_{1}}={x}_{2}$

${x}_{2}=\dot{{q}_{1}}\Rightarrow \dot{{x}_{2}}=-{M}^{-1}(1,:)\cdot C\cdot {\left[\begin{array}{cc}{x}_{2}& {x}_{4}\end{array}\right]}^{T}-{M}^{-1}(1,:)\cdot G+{M}^{-1}(1,:)\cdot \tau $${x}_{3}={q}_{2}\Rightarrow \dot{{x}_{3}}={x}_{4}$

${x}_{4}=\dot{{q}_{2}}\Rightarrow \dot{{x}_{4}}=-{M}^{-1}(2,:)\cdot C\cdot {\left[\begin{array}{cc}{x}_{2}& {x}_{4}\end{array}\right]}^{T}-{M}^{-1}(2,:)\cdot G+{M}^{-1}(2,:)\cdot \tau $

Suppose I would like to use a solver that takes as an argument the mass matrix (a MATLAB ODE solver in particular) and don't use its inverse because this will also simplify the computation of the jacobian (I intend to simulate a 7-DOF robotic manipulator after that so providing the mass matrix would be great). How can I write the initial equation (1) as a series of first order ordinary differential equations and be able to simulate its response by using some software solvers ?