I was wondering about what would be the most computationally efficient way in Julia to compute the zeros of the following equation involving Bessel functions derivatives:

J_m'(x)Y_m'(ax)-J_m'(ax)Y_m'(x)=0,

where a is a constant that can take values from (0,1] and m is an integer. I thought about using the asymptotic expressions to get an initial guess of the root and then improve it with a few Newton iterations, but I was wondering if there was a more efficient way to achieve this.

The SpecialFunctions package already has chainrules (e.g. Zygote) support for differentiating Bessel functions.

This is certainly the most common way to solve this kind of thing. If you want to find roots over and over (e.g. for different values of a), you can find roots at some set of points by Newton’s method and then interpolate. (Since it’s a smooth function, you can use exponentially accurate Chebyshev interpolation, e.g. via ApproxFun or FastChebInterp.)