I need an efficient way to numerically find the first $n$ positive roots $\lambda_n$ of the transcendental equation
$$ \dfrac{J_0 (\lambda_n r) Y_1 (\lambda_n) - J_1 (\lambda_n) Y_0 (\lambda_n r)}{J_0 (\lambda_n r) Y_0 (\lambda_n) - J_0 (\lambda_n) Y_0 (\lambda_n r)} = \dfrac{\delta}{\lambda_n}$$
where $\delta$ is a constant.
I have tried newtzero and chebfun in MATLAB. They are both taking longer than I'd like. The roots are needed for the solution of a PDE (transient heat transfer in one spatial dimension) with time varying boundary conditions using eigenfunction expansions. If anyone knows of a technique to determine the roots even a little faster, it would be much appreciated. Thank you.