Indeed, MassMatrix method isn't well documented, but this warning isn't hard to understand.
According to the tutorial Numerical Solution of Differential-Algebraic Equations, MassMatrix is a method that transforms the equation(s) to the form $M.x'=F(x,t)$, but what's undocumented is, the $M$ matrix must be a constant numerical matrix. What's constant numerical matrix? It's better to explain with a minimal example. Suppose we want to solve
$$t\ x'(t)-1=0,\ x(2)=3,\ t\in[2,4]$$
with MassMatrix method, how should we transform the ODE to the required $M.x'=F(x,t)$ form? Clearly the most straightforward transformation would be $(t).(x'(t))=(1)$ ({{t}}.{x'[t]} == {1}) i.e. $M=(t)$ (M = {{t}}). This is exactly what happens in NDSolve:
asso = NDSolve[{t x'[t] - 1 == 0, x[2] == 3}, x, {t, 2, 4},
Method -> EquationSimplification -> MassMatrix] // EvaluationData
NDSolve::ndmmc
asso["MessagesExpressions"][[1, 1, -1, 1]] // Normal
(* {{t}} *)
Remark
If I have to guess, NDSolve might have used something similar to
CoefficientArrays[{t x'[t] - 1 == 0}, {x'[t]}] // Last
to generate the mass matrix, but TraceInternal doesn't confirm it.
The $M$ involves symbolic variable $t$ (t) so it's not numerical. "But once the numeric ODE/DAE solver starts working, $t$ will be numerical at every time step!" Yeah, but it varies as time changes i.e. the matrix won't be constant. I'm not knowledge enough to explain why MassMatrix is designed not to allow such matrix, though.
"OK, then how to circumvent this? " In many cases, we can simply help the coefficient extractor a bit. For example, the example above can be transformed to
NDSolve[{x'[t] - 1/t == 0, x[2] == 3}, x, {t, 2, 4},
Method -> EquationSimplification -> MassMatrix]
"Wait, isn't this ridiculous? Why not directly set EquationSimplification -> Solve? That'll automatically transform the equation to x'[t] == 1/t!" Yeah, MassMatrix method doesn't show any advantage for this toy example. But when the system gets large (this happens quite a bit when symbolically discretizing PDE), the EquationSimplification -> Solve may be too expensive for symbolic pre-processing, while EquationSimplification -> Residual may be too weak to handle the resulting system, then MassMatrix method becomes an appreciable method. That's the reason I sometimes set this option when solving ODE systems generated by pdetoode. Typical example:
Why does NDSolve fail to solve the PDEs and spit out mconly warning?
Remark
For those who feel the toy system $(t).(x'(t))=(1)$ with only one
equation inside uncomfortabe, here's a slightly more complicated toy
example:
NDSolve[{{(t + 1) (x'[t] + y'[t]) == (t + 1) (Sin[t] - z'[t]),
y'[t] + z'[t] == Sin[t] + x[t],
x'[t] + z'[t] == y[t] + Cos[t]},
{x[0] == 1, y[0] == 1, z[0] == 1}},
{x, y, z}, {t, 0, 1},
Method -> {"EquationSimplification" -> "MassMatrix"}]
NDSolve::ndmmc
Notice this is essentially the documented example for MassMatrix, except that I've multiplied each side of the first equation
by (t + 1), yet this is enough to fool the pre-processor of
MassMatrix.
