I'm interested in a good reference/paper about how to handle *numerically* a mass-matrix system as 

\begin{align}
\mathbf{M}(t,y)\dot{y} =F(y,t)
\end{align}

I know that such a problem can be solved by using classical ODEsolvers in MatLab such as `ode15s`, and others. Of course, due to copiright issues, it's not possible to read the source code and, IMHO, there's a lack of reference.

I've found something interesting in Hairer & Wanner's books (both volumes), but honestly there's not so much stuff. 
In such problems, the non-singularity of $M(t,y)$ plays a fundamental role, and I'm quite interested in some techniques. 

- - - 
Suppose to have a *not-time-dependent* and *non-singular* mass matrix $M(t,y)=M(y)$ and a system 
\begin{align}
M(y)\dot{y} =F(y,t)
\end{align}

and say I want to compute the numerical solution with Backward Euler (or any implicit method, in order to avoid numerical instabilities)

So, one can write $\dot{y} = M^{-1}(y) f(y,t):= \tilde{F}(y,t)$ and formally we have

\begin{align}
y_{n+1}=y_n + \Delta t \tilde{F}(t_{n+1},y_{n+1})=y_n+ M(y_{n+1})^{-1} \Delta t F(t,y_{n+1})
\end{align}

and hence

\begin{align}
M(y_{n+1}) [y_{n+1}-y_n]=\Delta t F(t_{n+1},y_{n+1})
\end{align}

Another interesting problem is how to solve with Newton's nethod this non-linear system of equations. As far I know, a simplified Newton's method is used, but I can't found any reference about it.


I'd be very grateful if someone has references/hints or something else !