The answers to the question Bessel Functions Sturm-Liouville problem show that the eigenfunctions of the Helmholtz equation $$ (-\nabla^2+k^2)\psi(\mathbf r)=0 $$ in an annulus $a<r<b$, that have cylindrical separability of the form $\psi(r,\theta) = f(r)e^{im\theta}$ and which are subject to hard Dirichlet boundary conditions at both an inner and an outer ring, $$ f(a)=f(b)=0, $$ are linear combinations of Bessel functions of the first and second kinds, of the form $$ f(r) = N\bigg[Y_m(ka)J_m(kr)-J_m(ka)Y_m(kr)\bigg], $$ where the (square root of the) eigenvalue $k$ is fixed by requiring that the eigenfunction vanish at the outer ring: $$ f(b)/N= Y_m(ka)J_m(kb)-J_m(ka)Y_m(kb)=0 \tag{$*$} $$
I am looking for efficient and reliable ways to solve this equation, i.e. given $m$ and $b/a$, find the $k$ that corresponds to the eigenvalue.
For the circular problems, the zeroes of both $J_m(kr)$ and $Y_m(kr)$ are widely tabulated (example) and are often available as standard parts of scientific software packages (example), but the annular case isn't nearly as well-established. What are good ways to numerically solve that equation and reliably come up with the correct roots?