It's not a real good idea to go with an answer solely on the grounds that 'solve' can find one from an expression in a book. With books, there are definitely issues with the notation, who is using what, so you have to make sure you are solving the right equations. The E and H fields on each side of the boundary r = a are
E: G*J1(k1*a) = alpha*G*H1(k2*a)
H: G*(1/mu1)*k1*J0(k1*a) = alpha*G*(1/mu2)*k2*H0(k2*a)
where E is in the theta direction, around the cylinder, and H is in the z direction, along its length. (There is also a radial component of B that is matched automatically). G is an overall amplitude that doesn't matter for matching, alpha is an unknown constant to get the sizes to match up, and k1 and k2 are
k1 = sqrt(n1^2 k^2 - kz^2) k2 = sqrt(n2^2 k^2 - kz^2)
kz is your beta, which must be the same in both regions since the solution has to have the same z dependence down the length of the cylinder. k1 is your kappa, and k2 is almost your gamma. This is assuming that in your definition of k, k = 2pi / lambda, the wavelength is the free space version, lambda = c / f.
On the right, H1 and H0, which are really denoted H1(+) and H0(+) are certain Hanckel functions. (There is also an H1(-) and H0(-) ).
Dividing one equation by the other yields
mu1*J1(k1*a) / (k1*J0(k1*a)) = mu2*H1(k2*a) / (k2*H0(k2*a))
For a solution, kz has to be such that k1 is real and k2 is pure imaginary. Rewrite this as
k2 = i*g2 where g2 = +sqrt(kz^2 - n2^2 k^2)
and g2 is your gamma.
The fields in the outside region have to decay exponentially for large r, and the bessel functions K have that property. It turns out that
H0(i*g2*a) = -i*(2/pi)*K0(g2*a)
H1(i*g2*a) = -(2/pi)*K1(g2*a)
and after tossing those in we arrive at last at
mu1*J1(k1*a) / (k1*J0(k1*a)) = - mu2*K1(g2*a) / (g2*K0(g2*a))
Here the J's and K's are straight Matlab besselj and besselk.
Note the minus sign on the right. Various definitons of bessel functions include factors of (-)^n where n is the order of the bessel function. That messes with the minus signs.
This reply has gone on for a bit, but I beleive the process is absolutely necessary because of all the ways of defining bessel functions.
The larger the value of k*a, the more modes (roots) there will be. Your value of k*a is 40.54 which is large enough for several modes.
The code below shows where the roots are. I will leave it to you to use fzero or whatever method you choose to find roots.
a = 10e-6;
k = 2*pi/1.55e-6;
mu1 = 1;
n1 = 1.4;
mu2 = 1;
n2 = 1.33;
kzmax = n1*k;
kzmin = n2*k;
kz = linspace(kzmin,kzmax,100000);
kz(1) = ;
kz(end) = ;
k1 = sqrt(n1^2*k^2 - kz.^2);
g2 = sqrt(kz.^2 - n2^2*k^2);
C = mu1*besselj(1,k1*a)./(k1.*besselj(0,k1*a));
Cm = mean(abs(C));
C(abs(C)>3*Cm) = nan;
D = -mu2*besselk(1,g2*a)./(g2.*besselk(0,g2*a));
Dm = mean(abs(D));
D(abs(D)>3*Dm) = nan;
kzroot = 5.4617301149e6;
k1root = sqrt(n1^2*k^2 - kzroot^2);
k2root = sqrt(n2^2*k^2 - kzroot^2);
alpha = besselj(1,k1root*a)/besselh(1,1,k2root*a);
E1 = besselj(1, k1root*a)
E2 = alpha*besselh(1,1,k2root*a)
H1 = (1/mu1)*k1root*besselj(0, k1root*a)
H2 = alpha*(1/mu2)*k2root*besselh(0,1,k2root*a)