I want to solve a 'large' System of ODEs in the form:

$$\frac{d\mathbf{y}(t)}{dt}+M(t)\cdot\mathbf{y}(t)=0 \label{a} \tag{1}$$

Where M is an NxN Matrix, N is in the order of 600 and involves interpolating functions.

By just using the brute force approach:


where I just give Mathematica the raw list of flattened-out equations. Mathematica fails and suggests me to use:

Method -> {"EquationSimplification" -> "Residual"}

This leads to an overwhelming long initialization time before the algorithm then starts to solve the Equations (i tracked this by monitoring t). My Question are now:

  • How can I provide Mathematica directly a solvable equation-form without the need of the long initialization using Method -> {"EquationSimplification" -> "Residual"}. This Tutorial explains the method but i dont see how I can provide a form like $F(y',y,t)=0$ directly
  • What is the most efficient (or efficient) way to solve an equation like \ref{a}, where M involves interpolating function?
  • 1
    $\begingroup$ Should we assume that there are no symmetries or other additional structures in the equations or initial conditions? $\endgroup$ Oct 23, 2017 at 22:48
  • 1
    $\begingroup$ Have a look at in the section Transient PDEs. Even if this is for FEM it does deal with large sets of ODEs. $\endgroup$
    – user21
    Oct 23, 2017 at 23:36
  • $\begingroup$ Can M be diagonalized? $\endgroup$
    – bbgodfrey
    Oct 24, 2017 at 0:07
  • 1
    $\begingroup$ Have you tried the vector form of NDSolve? Something like NDSolve[{y'[t] + M[t] . y[t] ==0, y[0] == {..}}, y, {t, tstart, tend}]? $\endgroup$
    – Carl Woll
    Dec 26, 2017 at 21:37