Hi,

not sure whether this helps; it is certainly not the most elegant solution, but it can be generalised. I use a univariate ODE as a test case and I only perform the operation once. This integrates a simple exponential decay until we reach the value 0.1:

sol = NDSolve[{x'[t] == -2 x[t], x[0] == 1, WhenEvent[x[t] == 0.1, "StopIntegration"]}, x, {t, 0, 5}]

We then determine the domain of the interpolating function - otherwise it extrapolates when we plot it.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
domain = InterpolatingFunctionDomain[First[x /. sol]]

We can then plot the result and the threshold value:

Plot[{x[t] /. sol, 0.1}, {t, 0, domain[[1, 2]]}, PlotRange -> All]

which gives:

I hope that this helps to achieve what you want.

Cheers,

M.