Compare the two solution sets below: in the first one (`x == False`), there is a constant parameter 5 while in the second one (`x == True`), this parameter is introduced as a dependent algebraic variable `r[t] == 5` (converting the ODE to a DAE). A DAE often requires a different `NDSolve` method, here I used `"EquationSimplification" -> {Automatic, "SimplifySystem" -> True}` which messes up the solution order as is apparent from the plots. ClearAll[t, R, A, P, r, x]; eq = { R'[t] == -.1 R[t], A'[t] == .1 A[t], P'[t] == 0, If[x, r[t] == 5, Nothing], R[0] == 1, A[0] == 1, P[0] == 1 }; sol1 = Block[{x = False}, NDSolve[eq, {R, A, P}, {t, 0, 200}, DependentVariables -> {R, A, P}]]; sol2 = Block[{x = True}, NDSolve[eq, {R, A, P, r}, {t, 0, 200}, DependentVariables -> {R, A, P, r}, Method -> {"EquationSimplification" -> {Automatic, "SimplifySystem" -> True}}]]; Grid[{ LogPlot[#[t] /. sol1, {t, 0, 200}, PlotLabel -> #] & /@ {R, A, P}, LogPlot[#[t] /. sol2, {t, 0, 200}, PlotLabel -> #] & /@ {R, A, P, r} }] ![Mathematica graphics](https://i.sstatic.net/XjLp0.png) Can anyone reproduce this bug? Or is it supposed to be like this? (Mathematica 10.4.0.0, Windows 7, 64bits)