differential equations - When does NDSolve parallelize ODE system solving?

I've long believed that NDSolve cannot make use of multiple cores to solve ODE system, but things seem to be different at least since v12. Consider the following toy example:

$Version
(* "12.0.0 for Microsoft Windows (64-bit) (May 19, 2019)" *)

(* Tested on a 8-core machine. *)

eq = 2 y''[x] == y'[x] - 3 y[x] - 4;
ic = {y[0] == 5, y'[0] == 7};

With[{n = 4096}, sys = Table[{eq, ic} /. y -> y@i, {i, n}];
                 NDSolve[sys, y /@ Range@n, {x, 0, 10^3}]]

Seems that when n >= 4096, NDSolve automatically parallelizes. I'm pretty sure this isn't the case in v9.

So my question is:

1. 
   

   In which version was this optimization introduced?

   
   


2. 
   

   In what situation does the parallelization happen? More tests show the threshold depends on the specific system being solved. In certain situations (e.g. this) parallelization seems to never happen.

   
   


3. 
   

   Can we make NDSolve parallel for n < 4096? Is the parallelization controlled by any option?

   
   
   

Answer

The number 8192 is an autoparallelization threshold (see SystemOptions["ParallelOptions"]). There are 2*4096 = 8192 variables, y[k] and y'[k], in the system, perhaps a coincidence. The threshold can be lowered with SetSystemOptions, and one can see that when it is run, autoparallelization kicks in.

With[{ops = SystemOptions["ParallelOptions"]},
 Internal`WithLocalSettings[
  SetSystemOptions[
   "ParallelOptions" -> 
    "VectorParallelLengthThresholds" -> {8192, 8192, 8192, 8192, 
       8192}/2

   ],
  With[{n = 4096/2}, sys = Table[{eq, ic} /. y -> y@i, {i, n}];
   NDSolve[sys, y /@ Range@n, {x, 0, 10^3}]],
  SetSystemOptions[ops]
  ]]

This is a system arithmetic feature (vectorization), not an NDSolve feature per se.