I have been trying to solve numerically using MATLAB a system of ODEs, where when I introduce a state dependent function on the mass matrix (see below $T_{m1}(\theta)$), the residuals are suddenly increasing by a couple of orders of magnitude even when I use solvers such as ode23t which is designed to take strongly state dependent mass matrices.
To combat this I have chosen ode45 and made it run in a loop where the last point of the solution of the previous iteration acts as the initial value for the next. Also, for each iteration loops I set $T_{m1}(\theta)$ as a constant. This levels out the residuals but takes a tremendous amount of time for only 4 degrees of freedom (4x4 matrices that take about 30mins to solve). So I was wondering weather you'd have any advice of a numerical scheme that could be more efficient.
So here is the system of equations for the system:
$[\mathbf{J}-T_{m1}(\theta)]\ddot{\theta}+\mathbf{C}\dot{\theta}-T_{m2}(\theta)\dot{\theta}^2+\mathbf{T_l}\dot{\theta}^{1.7}+\mathbf{K}\theta=T_g(\theta)$
Where $\mathbf{J}$, $\mathbf{C}$, $\mathbf{K}$, $\mathbf{T_l}$, are matrices and $T_{m1}(\theta)$, $T_{m2}(\theta)$, $T_g(\theta)$, are vector valued functions (and obviously then $\theta$ is a vector)
Thanks for your help in advance!
KMT.
P.S.
What I mean by residuals is the result of the above equation when I shift $T_g(\theta)$ to the left hand side. In an ideal world it should be equal to zero, but as the numerical scheme accumulates errors, its not.
I have also tried all of MATLAB default solvers changing the tolerances as well but the residuals were still unacceptably high.
The functions of $\theta$ are all continuous sinusoidal functions with the same period.