Stopping NDSolve when encountering stiffness

Question

I am solving a differential equation for different initial conditions using ParametricNDSolveValue. I need to look at the value of the solution at some later point, and classify the nature of the solution. For certain starting values, the solution is well-behaved and oscillatory. For certain other starting values (which are smaller than the oscillatory ones), the solution rapidly rises to very large numbers, and Mathematica gives me a stiffness/singularity error:

ParametricNDSolveValue::ndsz: At r == 0.0004036933699289324`, step size is effectively zero; singularity or stiff system suspected. >>

I am wondering if there is a way to stop NDSolve the moment it encounters stiffness and somehow let me know that this happened. I essentially only need to know the value of the starting point at which it goes from oscillatory to stiff, with high precision. At such a value, the solution should end up at an unstable equilibrium point, without oscillations or rapidly rising to large numbers. Is there a way to detect stiffness (I am assuming the problem here is stiffness) without Mathematica trying to find an extrapolated solution?

Here's the equation I am working with:

Potential[x_] := 5.154462413581529*^7 - 256000000*(1 - x)^2 + 128000000*(1 - x)^4 - 1.030892482716306*^8*x^2 + 8.049495987048458*^8*(1 - x)^2*x^2 + 5.154462413581531*^7*x^4

test = ParametricNDSolveValue[{(Derivative[2][t][r] + (2/r)*Derivative[1][t][r] - D[Potential[x], x] /. x -> t[r]) == 0, t[10^(-12)] == d, 
 Derivative[1][t][10^(-12)] == 0}, t[0.5], {r, 10^(-12), 1}, {d}]; 

And this is the code I am using to check this in case anyone is interested:

C = 0.0005;
d = 0.002; 
Under = 1; 
Monitor[Quiet[While[C >= 10^(-18), While[Under == 1, d = d - C; If[Abs[test[d]] > 1 || d < 0, Under = 0]; ]; d = d + C; C = C/10; Under = 1; ]; ], d]

Answer 1

NDSolve stops itself when it believes that it has encountered stiffness. To take advantage of this, rewrite test as

test = ParametricNDSolveValue[{(Derivative[2][t][r] + (2/r)*
    Derivative[1][t][r] - D[Potential[x], x] /. x -> t[r]) == 0, 
    t[10^(-12)] == d, Derivative[1][t][10^(-12)] == 0}, t, {r, 10^(-12), 1}, {d}]

In other words, obtain the entire solution t, not just t[.5]. Then use, for instance,

test[.001]["Domain"][[1, 2]]
(* 0.000383406 *)

which shows that NDSolve detected stiffness at x == 0.000383406. On the other hand,

test[.002]["Domain"][[1, 2]]
(* 1. *)

shows that NDSolve did not encounter stiffness throughout the range of integration. Applying this to the code in the question that seeks the value of d for which stiffness first occurs gives

c = 0.0005; d = 0.002; Under = 1; Monitor[
Quiet[While[c >= 10^(-18), While[Under == 1, d = d - c; If[test[d]["Domain"][[1, 2]] < 1, 
    Under = 0];]; d = d + c; c = c/10; s = d; Under = 1;];], d]
NumberForm[s, 16]
(* 0.001701281449747991 *)

Note that c is used here instead of C, which is a reserved symbol in Mathematica.

Addendum: Faster Solution

The solution above, patterned after that in the question, requires 108 sec on my PC to obtain an answer. The following, which brackets the solution and successively reduces the bounds by a factor of two on each iteration,

dl = c; du = d; 
Do[dt = (dl + du)/2; If[Quiet@test[dt]["Domain"][[1, 2]] < 1, dl = dt, du = dt], {i, 20}]
NumberForm[du, 16]
(* 0.001701281449747996 *)

takes 0.000112592 sec.