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:


(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:

Solution to h[r]

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
  • $\begingroup$ What is Potential[a, b]? $\endgroup$
    – Sumit
    May 20, 2016 at 19:24
  • $\begingroup$ I have edited the question to include this. $\endgroup$
    – Gowri
    May 20, 2016 at 19:38
  • $\begingroup$ 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 $\endgroup$
    – george2079
    May 20, 2016 at 20:51
  • $\begingroup$ 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 ? $\endgroup$
    – Gowri
    May 20, 2016 at 20:51

2 Answers 2


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

enter image description here

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}]

enter image description here

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.

enter image description here

  • $\begingroup$ +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. $\endgroup$
    – Sumit
    May 20, 2016 at 21:23
  • $\begingroup$ Thanks for the answer! Do you have any idea why it can't solve past that point? $\endgroup$
    – Gowri
    May 21, 2016 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}]

enter image description here

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


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}

  • $\begingroup$ 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. $\endgroup$
    – Gowri
    May 20, 2016 at 20:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.