'Returns the system\'s equations of motion in first order form. The
output is the right hand side of::
x\' = |q\'| =: f(q, u, r, p, t)
|u\'|
The right hand side is what is needed by most numerical ODE
integrators.
Parameters
inv_method : str
The specific sympy inverse matrix calculation method to use. For a
list of valid methods, see
:meth:`~sympy.matrices.matrices.MatrixBase.inv`'

def rhs(self, inv_method=None):

    rhs = zeros((len(self.q) + len(self.u)), 1)
    kdes = self.kindiffdict()
    for (i, q_i) in enumerate(self.q):
        rhs[i] = kdes[q_i.diff()]
    if (inv_method is None):
        rhs[len(self.q):, 0] = self.mass_matrix.LUsolve(self.forcing)
    else:
        rhs[len(self.q):, 0] = (self.mass_matrix.inv(inv_method, try_block_diag=True) * self.forcing)
    return rhs
