I am trying to use a solution given by Michael E2 in this topic:

Zeros[z_] = (BesselJ[0, z] - BesselJ[2, z])*(3*BesselI[1, z] + 
  BesselI[3, z]) - (BesselI[0, z] + 
  BesselI[2, z])*(3*BesselJ[1, z] - BesselJ[3, z]);

(zeros = First @ Last @ Reap @ Quiet @ NDSolve[
   {y'[x] == Zeros'[x], y[0] == Zeros[0], 
    WhenEvent[y[x] == 0, Sow[FindRoot[Zeros[z], {z, x}]]]}, 
   y, {x, 0, 40}, AccuracyGoal -> 1, PrecisionGoal -> 1]) // AbsoluteTiming

-- to find roots for my equation:

$$ J_1(a/z)\, Y_0(b/z)-J_0(a/z)\,Y_1(b/z)=0, \quad a=10, \;b=1, $$

but get an error

First::first: "{} has a length of zero and no first element."

The problem is that "y(x)=0" event doesn't happen. I understand that the functions are quite different, so I perform the following:

$$ J_1(a/z)\; Y_0(b/z)-J_0(a/z)\; Y_1(b/z)=0 \rightarrow\\ \rightarrow \left( J_1(a\tilde z)\; Y_0(b\tilde z)-J_0(a\tilde z)\; Y_1(b\tilde z) \right)\; e^{-100/\tilde z}=0, ~\tilde z = 1/z $$

to make my function to be more like the original one, but still get the error. Could anyone help, please?

My code (I mean, the original code with my function) is below.

Sigma = 1;
xi1 = 10;
Zeros[z_] = 
  (BesselJ[1, xi1*z]*BesselY[0, Sigma*z] - 
   BesselJ[0, Sigma*z]*BesselY[1, xi1*z])*Exp[-100/z];
(zeros = First@
  Quiet@NDSolve[{y'[x] == Zeros'[x], y[0] == Zeros[0], 
 WhenEvent[y[x] == 0, Sow[FindRoot[Zeros[z], {z, x}]]]}, 
y, {x, 0, 10}, AccuracyGoal -> 10, PrecisionGoal -> 10]) // AbsoluteTiming
  • $\begingroup$ The first code you posted works OK here. Try restarting Mathematica. If it doesn't work, post the version of Mathematica you're using $\endgroup$ Nov 20, 2013 at 20:03
  • $\begingroup$ @belisarius First code works good. The problem is that I can't get it working with another function, even if I make that function look like the first one (more or less). See second code. My version of Mathematica is, 64-bit. $\endgroup$
    – ApplMath
    Nov 23, 2013 at 10:51

1 Answer 1


Solved it! Here's the code.

xi1 = 10;
Sigma = 1;
f3[z_] = BesselJ[1, xi1/z]*BesselY[0, Sigma/z] - 
   BesselJ[0, Sigma/z]*BesselY[1, xi1/z];

Timing[roots3 = 
 NDSolve[{x'[z] == f3'[z], x[1] == f3[1]}, {x}, {z, 15, 0.1}, 
  Method -> {"EventLocator", "Event" -> x[z], 
    "EventAction" :> Sow[z]}, MaxStepSize -> 0.001, 
  MaxSteps -> 1000000]]][[2, 1]]]


f3 /@ roots3

shows that the roots are quite accurate. Sure it doesn't find all the roots in the interval and skips some of them, but one can easily increase the accuracy.

Note: I used my initial function.


Your Answer

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

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