Solving one unknown variable in the matrix with matlab symbol operation

clc
clear
epsilon_Ag = (-11.046476 + 4.43688i);
epsilon_vac = 1;
c_light =3 * 10^8;
omega =2*pi*5.765e14;%wavelength = 520nm
k1 = omega / c_light*sqrt(epsilon_vac);
k2 = omega / c_light*sqrt(epsilon_Ag);%the wave vector of Ag
u_0 =1;%or 4*pi*10^(-7);
r =5e-8;%50nm
syms k_p;
k1o =sqrt(k1^2-k_p^2);
k2o =sqrt(k2^2-k_p^2);
besselh0=besselj(0,k1o*r)+1i*bessely(0,k1o*r);%hankel
diff_h0 =-besselj(1,k1o*r) - bessely(1,k1o*r)*1i;

M =[k1o^2/k1^2*besselh0 -k2o^2/k2^2*besselj(0,k2o*r);1i/omega/u_0*k1o*diff_h0 1i/omega/u_0*k2o*besselj(1,k2o*r)];%matrix of M
k_para0 =vpasolve(det(M),k_p);%solve det M =0
k_para0 =double(k_para0);

k1o_0 =sqrt(k1^2-k_para0^2);
k2o_0 =sqrt(k2^2-k_para0^2);

I want to solve the det M=0 to get the value of k_p, but the value of k_p i got is larger than k1(about 10^6 times). In fact,the value of k_p is smaller than k1, or we can say that they are of the same order of magnitude. Can u give me some suggestions?

[guodong zhu]

[guodong zhu]

I try to solve the determinant with mathematica, the computation consume a lot of running memory and I haven't got the result yet.