# Difficulty finding roots of an interpolation function from NDSolve

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

• What is Potential[a, b]? – Sumit May 20 '16 at 19:24
• I have edited the question to include this. – Gowri May 20 '16 at 19:38
• Be sure to look at the actual domain of the resulting solution (Look at the interpolation function display). When I run this it halts at 0.0472 with a step size zero error. If you ignore that error and plot the result past r=0.0472 you are not seeing a solution at all but just what results from extrapolating the interpolation function. h is strictly positive and blows up at 0.0472 so it never gets to your WhenEvent – george2079 May 20 '16 at 20:51
• I am not sure I understand what you said about the domain. I am using NDSolve in the range {r,10^-18,1}. Why would it not calculate beyond r=0.0472 ? – Gowri May 20 '16 at 20:51

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.

• +1 george. I put a quiet when I saw 10^-18 in the main function. I didn't understand that the extrapolation can be a rubbish. – Sumit May 20 '16 at 21:23
• Thanks for the answer! Do you have any idea why it can't solve past that point? – Gowri May 21 '16 at 15:16

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}

• I get that I can find multiple roots by changing the initial point, but I was looking for a way to get all the roots in a particular range without ambiguity. Although the FindRoot method might work in my particular case, I was curious if there was a more efficient and better way of doing this. – Gowri May 20 '16 at 20:45