(* 数値データの準備 *)

Clear[bj,dbj,bi,dbi,mu,a,b,c,d,sol,kmmin,kmn,lambda]

bj[m_,x_]:=BesselJ[m,x];

dbj[m_,x_]:=(BesselJ[-1+m,x]-BesselJ[1+m,x])/2;

bi[m_,x_]:=BesselI[m,x];

dbi[m_,x_]:=(BesselI[-1+m,x]+BesselI[1+m,x])/2;

mu=N[3/10,100];

a[m_,k_]:=(mu-1) (k dbj[m,k]-m^2 bj[m,k])-k^2 bj[m,k];

b[m_,k_]:=(mu-1) (k dbi[m,k]-m^2 bi[m,k])+k^2 bi[m,k];

c[m_,k_]:=m^2 (mu-1) (k dbj[m,k]-bj[m,k])-k^3 dbj[m,k];

d[m_,k_]:=m^2 (mu-1) (k dbi[m,k]-bi[m,k])+k^3 dbi[m,k];

sol[m_,near_]:=k/.FindRoot[a[m,k]/b[m,k]==c[m,k]/d[m,k],{k,near},WorkingPrecision->20];

kmmin[0]=sol[0,3.0];

kmmin[1]=sol[1,4.5];

kmmin[2]=sol[2,2.5];

kmmin[m_]:=kmmin[m]=sol[m,kmmin[m-1]+1];

Table[kmmin[m],{m,0,50}];

kmn[0,0]=0;

kmn[1,0]=0;

kmn[0,1]=kmmin[0];

kmn[1,1]=kmmin[1];

Table[kmn[m,0]=kmmin[m],{m,2,50}];

kmn[m_,n_]:=kmn[m,n]=sol[m,kmn[m,n-1]+Pi]

kmns=Table[kmn[m,n],{m,0,50},{n,0,20}];

Export["kmns.dat",kmns];

lambdamn[m_,n_]:=lambdamn[m,n]=If[kmn[m,n]>0,-a[m,kmn[m,n]]/b[m,kmn[m,n]],0]

lambdamns=Table[lambdamn[m,n],{m,0,50},{n,0,20}];

Export["lambdamns.dat",lambdamns];