Labeling solutions of an Eigenvalue equation involving Bessel functions

I've worked with a very similar equation quantization condition. Basically this is a root-finding problem, but since there are many roots and you want exactly the $n^\textrm{th}$, FindRoot is not very helpful. The approach which worked for me is to use NDSolve along with Sow and Reap. This question spells out the basics.

I'd try something like this:

quantcond[m_Integer, k_, R1_, R2_] := BesselJ[m, k R1] BesselY[m, k R2] - BesselJ[m, k R2] BesselY[m, k R1]

quantcondZero[m_?NumericQ, R1_?NumericQ, R2_?NumericQ, n_Integer, opts : OptionsPattern@NDSolve] := Module[
  {r = 0,
   min = 1/(1000 Min[R1, R2]),
   max = n 1000/Min[R1, R2]},

  Reap[
    NDSolve[
     {f'[k] == D[quantcond[m, k, R1, R2], k],
      f[min] == quantcond[m, min, R1, R2],
      WhenEvent[f[k] == 0, Sow[k]; r++],
      WhenEvent[r > n, "StopIntegration"]},
     f,
     {k, min, max},
     opts
    ]
   ][[2, 1]]
]

The function quantcondZero will return the first n zeros of the function quantcond.

Example

quantcondZero[1, 20^-3, 30^-3, 9]

{3160.94, 6293.06, 9431.39, 12571.3, 15711.9, 18852.9, 21994., 25135.2, 28276.5}

The roots look almost like multiples of π, which is the expected behavior.

Notes

1. The main difficulty lies in scaling. Here, assuming that $R_1$ and $R_2$ are roughly the same scale, your zeros appear to scale roughly inversely with the size of $R_1$ and $R_2$, hence the setting for min and max, the limits of the integration. Depending on the behavior of the function, these may need to be adjusted to ensure you capture all of the zeros.

I typically solve these problems by harnessing the power of Plot. Taking

fn[m_, k_, R1_, R2_] := BesselJ[m, k R1] BesselY[m, k R2] - BesselY[m, k R1] BesselJ[m, k R2]

we first Plot the function as

p1 = Plot[fn[1, k, 1, 2], {k, 0, 50}
  , PlotRange -> All
  , MeshFunctions -> {#2 &}
  , MeshStyle -> {Red, PointSize[0.02]}
  , Mesh -> {{0.}}]

yielding

We extract the roots and refine them using FindRoot:

pts = Sort@Cases[Normal@p1, Point[a_] :> First@a, Infinity]
k /. FindRoot[fn[1, k, 1, 2] == 0, {k, #}] & /@ pts

We can wrap this all up into a Module with the following code (borrowing Virgil's name from their answer):

quantcondZero[m_?NumericQ, R1_?NumericQ, R2_?NumericQ, n_Integer, x_?Positive] := Module[
  {plt, pts}
  , plt = Plot[BesselJ[m, k R1] BesselY[m, k R2] - BesselY[m, k R1] BesselJ[m, k R2]
      , {k, 0, x}, PlotRange -> All
      , MeshFunctions -> {#2 &}, Mesh -> {{0.}}]
  ; pts = Sort@Cases[Normal@plt, Point[a_] :> First@a, Infinity]
  ; k /. FindRoot[eqn[1, k, 1, 2] == 0, {k, #}] & /@ pts[[1 ;; Min[n, Length@pts]]]
  ]

Because one can't tell a priori how far out to plot the function, I have to include it as an argument to the function, and it makes it not completely automated. So, for instance, if we want the first 5 roots with $m=2$, we can start by doing

quantcondZero[2, 1, 2, 5, 10]
(* {3.40692, 6.42777, 9.52285} *)

It doesn't give us enough, so we try a larger x:

quantcondZero[2, 1, 2, 5, 20]
(* {3.40692, 6.42777, 9.52285, 12.6404, 15.7673} *)

---

Important note: this method is not the most robust. If you need a lot of zeroes, it is likely that the zeroes will get so close together that the Plot might miss them. You will need to do things like manually increase the number of PlotPoints in order to get it to work.

My goal was to give an answer using FindRoot. In addition, I want to be able to uniquely identify the roots by the node number. I believe this is an important part of the original question, because it refers specifically to "labeling" the solutions. The method below achieves this unique association of the label n with the number of zeros. For a given n, it does not require finding all roots up to n and counting them. This is where it differs fundamentally from the other answers.

The main issue is to find the starting points for FindRoot reliably without duplicating or skipping roots. For large enough m, this can be done very nicely by relying on BesselJZero, but this doesn't work for all roots.

Therefore, I decided to modify my answer and not rely on BesselJZero at all. The alternative is to use an asymptotic approach based on the WKB approximation. By doing this approximation with a rescaled coordinate $\xi \equiv \log(\rho)$, you get a very accurate formula for the starting points that works for large and small m. I compared it to the solution available in the package linked in the question, and found it promising that my function doesn't duplicate roots whereas the library function does so:

BesselJYJYZeros[20, .5, 10]

(*
==> {25.4172, 25.4172, 29.9645, 34.0415, 38.1114, 42.5018, \
47.2888, 52.4004, 57.7525, 63.2819}
*)

You can see that the first root is duplicated, showing that the search was started with bad initial conditions.

Here is my approach:

Clear[annularZero, action, sWKB];

action[m_, k_, ρa_] := 
 Sqrt[k^2 - m^2] - Sqrt[-m^2 + k^2 ρa^2] + 
  m ArcTan[(m (-Sqrt[k^2 - m^2] + Sqrt[-m^2 + k^2 ρa^2]))/(
    m^2 + Sqrt[k^2 - m^2] Sqrt[-m^2 + k^2 ρa^2])]

sWKB[m_, k_?NumericQ, r_] :=
  Re[action[m, k, Max[m/k, r]]];

Options[annularZero] = {"ShowPlot" -> False, 
   WorkingPrecision -> $MachinePrecision};

annularZero[m_, ρa_, n_, OptionsPattern[]] := Module[
   {k, det, start1, start2, refine = False, firstPass, prec, showPlot},
   (* For large m it helps to increase WorkingPrecision if FindRoot \
complains: *)
   prec = OptionValue[WorkingPrecision];
   (* With showPlot True, display the function whose roots we seek: *)
      showPlot = OptionValue["ShowPlot"] =!= False;
   (* Use WKB approximation to find two starting values for exact \
root search. No need for perfection, 
   so suppress warnings: *)
   {start1, start2} = 
    Table[
     Chop[k /.
       Quiet@FindRoot[
         sWKB[m, k, ρa] - Pi (n + δ), {k, 
          m/ρa}]], {δ, {-1/4, 0}}];
   (* Want the n-th exact root of det: *)
   det = SetPrecision[
     BesselJ[m, k ρa] BesselY[m, k] - 
      BesselJ[m, k] BesselY[m, k ρa], prec];
   (* Since det has infinitely many roots, 
   optionally show plot of interval around desired region: *)
   If[showPlot,
    Print[
     Plot[det, {k, Max[0, start1 - 10], start2 + 10}, 
      Epilog -> {Red, PointSize[.02], 
        Point[{{start1, 0}, {start2, 0}}]}, PlotRange -> {-1, 1}]]];
   (* FindRoot with two starting values to bracket the desired root:*)
      Quiet[
    Check[
     firstPass = Chop[k /. First@FindRoot[det, {k, start1, start2}]], 
     refine = True]];
   (* If previous search issues warnings, 
   try to improve last result using Automatic method: *)
   If [refine, If[showPlot, Print["refining zero"]]; 
    Chop@N[k /. 
       First@FindRoot[det, {k, firstPass}, WorkingPrecision -> prec]],
     firstPass]
   ];

For now, I simplified life for myself by requiring that the second argument ρa must be less than 1. This can always be achieved by rescaling the original problem to appropriate length units. I may drop this restriction once the reliability has been tested some more. Here is the comparison to the test case above:

Table[annularZero[20, .5, i], {i, 10}]

(*
==> {25.4172, 29.9645, 34.0415, 38.1114, 42.5018, 47.2888, \
52.4004, 57.7525, 63.2819, 68.9443}
*)

Instead of duplicating the first root, I now get a new additional root at the end of the list. The other results are in complete agreement with the old add-on package.