For j = 1, figure 1 shows the general behavior of the function. Since it's a plot of log(abs(fun)) there can be no zero crossings but you can see the attempts to head for zero. There are a lot of roots but no roots below mu = 1 which is also true for j = 2,3,4.
Since it's an analytic function (no noise) and the roots are very close to equally spaced (which will become evident in figure 2), a fairly simple method can be used to find the approximate locations of the zeros: convert the function to a square wave and find the discontinuities.
With fzero. rather than feed it an approximate location for a zero, it is safer to give it an upper and lower bound where there is sure to be exactly one zero crossing in between. The code finds the points halfway between the approximate zeros, and adds 1 at the beginning since there is no zero to the left of that.
Plot 2 shows a result for j = 1. The function is dropping in amplitude so rapidly in that region that it helps to multiply by a factor of mu.^15 just to make it easy to see. That of course does not affect the location of the zeros.
Results for j = 2,3,4 are similar.
R_i = 3;
R_o = 13;
modes = [19,18,17,16,15,14];
j = 1;
m = modes(j);
k = modes(j+2);
fun = @(mu) (0.5*(besselj(m,mu*R_o)-besselj(k,mu*R_o))).*(0.5*(bessely(m,mu*R_i)-bessely(k,mu*R_i))) ...
mu = .01:.01:50;
fin = find(abs(diff(sign(fun(mu))))>1);
mu0approx = mu(fin);
mu1 = [1 (mu0approx(1:end-1)+mu0approx(2:end))/2];
mu0 = ;
for n = 1:length(mu1)-1
mu0 = [mu0 fzero(fun,mu1(n:n+1))];
mu = 0:.001:6;
mu06 = mu0(mu0<6);