numerical integration - Why does NDSolve need to solve for the derivatives if the equations are already explicitly solved?

Consider the following test case:

mMax = 50;
Ffuncs[r_] = Table[F[n, m][r], {n, mMax}, {m, mMax}];
NDSolve[Flatten@{
    Thread[Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]],
    Thread[Flatten@Ffuncs[0] == Table[0, {mMax*mMax}]]

    }, Flatten@Ffuncs[r], {r, 0, 2}];

When I set mMax to 49 or less, it works without any problems. But with 50 or higher, I get

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

That was Mathematica 11. Mathematica 9 gives a better message and a better suggestion, which actually works, unlike the suggestion above, which takes forever to complete the solution and eventually crashes the kernel. Here's Mathematica 9's message:

NDSolve::ndsdtc: The time constraint of 1.` seconds was exceeded trying to solve for derivatives, so the system will be treated as a system of differential-algebraic equations. You can use Method->{"EquationSimplification"->"Solve"} to have the system solved as ordinary differential equations. >>

But anyway, I've already gave NDSolve a system of equations explicitly solved for derivatives! Why does it still need to solve for derivatives?

Answer

In V10 and higher, you can increase the time limit on Solve with the system option "NDSolveOptions" -> "DefaultSolveTimeConstraint" -> time:

With[{opts = SystemOptions[]},
 Internal`WithLocalSettings[
  SetSystemOptions["NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 10.`],
  sol = NDSolve[
     Flatten@{Thread[
        Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]], 
       Thread[Flatten@Ffuncs[0] == Table[0, {mMax*mMax}]]}, 

     Flatten@Ffuncs[r], {r, 0, 2}];,
  SetSystemOptions[opts]
  ]]
(* Spurious message: SystemOptions::noset: ..SystemOption PreemptiveCheckUseThreads..  *)

sol[[1, 1]]
(*  F[1, 1][r] -> InterpolatingFunction[{{0., 2.}}, <>][r]  *)

I'm not sure what takes NDSolve[] so long, but it's much faster than Solve[]:

Solve[Thread[Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]], 

   Flatten[D[Ffuncs[r], r]]]; // AbsoluteTiming
(*  {51.1422, Null}  *)