I'm trying to reset the value of functions $u_{i}$ when they reach a threshold $θ$. Below the threshold, they evolve according to a simple differential equation. The $r_{i}$ resets the derivative to zero for a while (refractory period of a neuron), but I doubt that's relevant to my question.

NDSolve is really slow when the number of events is increased (which can be done by decreasing recovery). I don't understand why an artificial discontinuity would make this code slower. Here is a minimal example where all $u_{i}$ are the same:

recovery = 0.01; θ = 0.5; n = 50; T = 2;
s = ConstantArray[{}, n];

AbsoluteTiming[
 sol = NDSolve[Table[With[{i = i}, {
  Subscript[u, i]'[t] == (-Subscript[u, i][t] + 10^3) Boole[Subscript[r, i][t] < t], 
  Subscript[r, i][0] == 0, 
  Subscript[u, i][0] == 0, 

  WhenEvent[
   Subscript[u, i][t] > θ, {Subscript[u, i][t] -> 0, 
   AppendTo[s[[i]], t], 
   Subscript[r, i][t] -> t + recovery}]}],
 {i, n}], 

  Table[{Subscript[u, i], Subscript[r, i]}, {i, n}] // Flatten, {t, 0, T}, 
  DiscreteVariables -> Table[Subscript[r, i], {i, n}]];
 ]

Length[s // First] (*number of events*)

(*{65.485, Null}*)
(*191*)

What could be causing this, and is there any way to make it faster?