Difficulty finding roots of an interpolation function from NDSolve

Question

I have been trying to find the point at which one of the solutions of a system of two ODEs crosses zero. I used the method suggested in this answer to a previous question, which seemed to be the most efficient:

https://mathematica.stackexchange.com/a/16448/38989

(PS- I had tried FindRoot, but it was not useful as it finds only one root).

My problem is that the solution I am looking at crosses zero and falls very rapidly to large negative numbers. The above code doesn't seem to work in this case. Here is the code I used:

DiffEq1 = Derivative[2][h][r] + (2/r)*Derivative[1][h][r] - D[Potential[a, b], a] /. {a -> h[r], b -> s[r]}; 
DiffEq2 = Derivative[2][s][r] + (2/r)*Derivative[1][s][r] - D[Potential[a, b], b] /. {a -> h[r], b -> s[r]}; 

Reap[NDSolve[{DiffEq1 == 0, DiffEq2 == 0, h[10^(-18)] == 225, s[10^(-18)] == 10, 
    Derivative[1][h][10^(-18)] == 0, Derivative[1][s][10^(-18)] == 0, 
    WhenEvent[h[r] == 0, Sow[s[r]]]}, {h, s}, {r, 10^(-18), 1}]]

The differential equations have a $2/r$ factor, which is why I am starting from $10^{-18}$. (This singularity is not an actual problem since the initial first derivative is 0).

When I run this, I don't get any solutions at all, even though it does work for some other initial starting points. I know that the solution to h[r] does cross zero at least once, as it starts from a positive value, and it ends up at a very large negative value:

I am struggling to find the reason for this. Does the Event Locator method have a problem with solutions that blow up and go to large numbers?

Edit: Potential[a,b] is a degree 4 polynomial in a and b. If I put in the numbers, this is what I get:

Potential[h_,s_]:= 1.1868147535541435*^8 - (15625*h^2)/4 + (15625*h^4)/524288 - 12768*s^2 + (h^2*s^2)/2 + 0.34340200000000004*s^4

Answer 1

I guess this is worth explaining a little better. Here is what I get running the exact code:

You can see the error message and the domain say the integration stopped around 0.0473.

The trouble arises here, if you don't notice that and plot over your specified range

 Plot[h[r] /. First@sol, {r, 10^-18, 1}]

you get a plot, but its just garbage past the error point. Why this is allowed baffles me. Note if you request a specific value, eg h[.05] /. First@sol you get a warning:

"Input value {0.05`} lies outside the range of data in the interpolating function. Extrapolation will be used. "

For some reason Plot does not give that warning.

Here is the actual solution:

 Plot[h[r] /. First@sol, {r, 
    10^-18, (h /. First@sol)["Domain"][[1, 2]]}]

which you see is always positive.

Answer 2

This is what I get with MMA10 (Linux).

{solh, sols} = {h, s} /. 
 NDSolve[{DiffEq1 == 0, DiffEq2 == 0, h[10^(-18)] == 225, 
 s[10^(-18)] == 10, Derivative[1][h][10^(-18)] == 0, 
 Derivative[1][s][10^(-18)] == 0, 
 WhenEvent[h[r] == 0, Sow[s[r]]]}, {h, s}, {r, 10^(-18), 1}][[1]];

Plot[{solh[x], sols[x]}, {x, 0, .35}]

x /. FindRoot[solh[x] == 0, {x, 0.3}]

0.0472792

You can find more zeros by changing the initial guess.

Table[x /. FindRoot[solh[x] == 0, {x, x0}], {x0, 0, 0.1, 0.01}]

{0., 0.0355549, 0.0366415, 0.03612, 0.0380459, 0.0472792, 0.0472792, \ 0.0472792, 0.0472792, 0.0472792, 0.0472792}