SetDirectory[NotebookDirectory[]];
<<MaTeX`
texStyle={};
SetOptions[MaTeX,"BasePreamble"->{"\\usepackage{amsmath}","\\usepackage{xcolor}","\\usepackage{fourier}","\\usepackage{ebgaramond}"},FontSize->11];

frame[legend_]:=Framed[legend,FrameStyle->Thin,RoundingRadius->10,FrameMargins->0];

Clear[findAllRoots]
(* from https://mathematica.stackexchange.com/questions/16439/find-all-roots-of-an-interpolating-function-solution-to-a-differential-equation/16444#16444 *)
SyntaxInformation[findAllRoots]={"LocalVariables"->{"Plot",{2,2}},"ArgumentsPattern"->{_,_,OptionsPattern[]}};
SetAttributes[findAllRoots,HoldAll];

Options[findAllRoots]=Join[{"ShowPlot"->False,PlotRange->All},FilterRules[Options[Plot],Except[PlotRange]]];

findAllRoots[fn_,{l_,lmin_,lmax_},opts:OptionsPattern[]]:=Module[{pl,p,x,localFunction,brackets},localFunction=ReleaseHold[Hold[fn]/. HoldPattern[l]:>x];
If[lmin≠lmax,pl=Plot[localFunction,{x,lmin,lmax},Evaluate@FilterRules[Join[{opts},Options[findAllRoots]],Options[Plot]]];
p=Cases[pl,Line[{x__}]:>x,Infinity];
If[OptionValue["ShowPlot"],Print[Show[pl,PlotLabel->"Finding roots for this function",ImageSize->200,BaseStyle->{FontSize->8}]]],p={}];
brackets=Map[First,Select[(*This Split trick pretends that two points on the curve are "equal" if the function values have _opposite _ sign.Pairs of such sign-changes form the brackets for the subsequent FindRoot*)Split[p,Sign[Last[#2]]==-Sign[Last[#1]]&],Length[#1]==2&],{2}];
x/. Apply[FindRoot[localFunction==0,{x,##1}]&,brackets,{1}]/. x->{}]

Coefficients of Weyl’s expansion in various settings

αλ[λ_,μ_]:=μ/(λ+2μ);
λα[α_,μ_]:=μ(1/α-2);
CWeyl[d_,λ_, μ_]:=((λ+2 μ)^(-d/2)+(d-1) μ^(-d/2)) /(4  Pi)^(d/2)/Gamma[1+d/2];
BDirLiu[d_,λ_, μ_]:=-((d-1)/μ^((d-1)/2)+1/(λ+2μ)^((d-1)/2))/4/(4 Pi)^((d-1)/2)/Gamma[1+(d-1)/2];
BDir[d_,α_, μ_]:=-μ^((1-d)/2)/(2^(d+1)Pi^((d-1)/2)Gamma[(1+d)/2]) (4(d-1)/Pi NIntegrate[τ^(d-2)ArcTan[Sqrt[(1-α τ^(-2))(τ^(-2)-1)]],{τ,Sqrt[α],1}]+α^((d-1)/2)+d-1);
BDirSaVa2[α_, μ_]:=1/(4Pi Sqrt[μ])(-1-Sqrt[α]-4/Pi NIntegrate[ArcTan[Sqrt[(1-α/ζ^2)(1/ζ^2-1)]],{ζ, Sqrt[α],1}]);
BDirSaVa3[λ_, μ_]:=-1/(16Pi)(3λ^2+13λ μ +16μ^2)/(λ^2 μ+5λ μ^2 +6μ^3);
γR[α_]:=Sqrt[Min[x/.Solve[(x^3)-8 (x^2)+16 (-1+α)+8(3-2 α) x ==0,x, Reals]]];
Bfree[d_,α_, μ_]:= μ^((1-d)/2)/(2^(d+1)Pi^((d-1)/2)Gamma[(1+d)/2]) (4(d-1)/Pi NIntegrate[τ^(d-2)ArcTan[(τ^(-2)-2)^2/(4Sqrt[(1-α τ^(-2))(τ^(-2)-1)])],{τ,Sqrt[α],1}]+α^((d-1)/2)+d-5+4γR[α]^(1-d));
BfreeSaVa2[α_, μ_]:=1/(4Pi Sqrt[μ])(4/γR[α]-3+Sqrt[α]+4/Pi NIntegrate[ArcTan[(2-1/ζ^2)^2/(4Sqrt[(1-α/ζ^2)(1/ζ^2-1)])],{ζ, Sqrt[α],1}]);
BfreeSaVa3[λ_, μ_]:=1/(16Pi) (3(λ+2μ)^2-3(λ+2μ)μ+2μ^2)/((λ+2μ)μ(λ+μ));

Lame operator, L, is the operator of plane elasticity written in polar coordinates. For vectors, the fist component is along the radius.

Lame = Function[u, (μ Curl[Curl[Flatten[{u,0}],{r,ϕ,z},"Cylindrical"],{r,ϕ,z},"Cylindrical"]-(λ+2μ) Grad[Div[Flatten[{u,0}],{r,ϕ,z},"Cylindrical"],{r,ϕ,z},"Cylindrical"])[[1;;2]]];
Lame[{(Subscript[u, 1])[r,ϕ],(Subscript[u, 2])[r,ϕ]}]

The boundary conditions for the free boundary and for the  Dirichlet problems are set up with

T[u_]:=({λ Div[u, {r, ϕ}, "Polar"],0}+2μ D[u,r]+μ({1,0,0}\[Cross]Curl[Flatten[{u,0}],{r,ϕ,z},"Cylindrical"])[[1;;2]]/.{r->1})//Simplify;
TDir[u_]:={u[[1]],u[[2]]};

T[{(Subscript[u, 1])[r,ϕ],(Subscript[u, 2])[r,ϕ]}]
TDir[{(Subscript[u, 1])[r,ϕ],(Subscript[u, 2])[r,ϕ]}]

Solutions of Lame system for Λ=0: there are 4 solutions, {r,0}, {0,r}, {1, I} Exp[I ϕ], {1, -I} Exp[-I ϕ] 
This is irrelevant for the discussion and can be skipped

Collect[Lame[{a r,c r}], r, Simplify]
T[{a r,c r}]

Collect[Exp[-I ϕ]Lame[{b,I b} Exp[I ϕ]], r, Simplify]
T[{b,I b} Exp[I ϕ]]

Collect[Exp[I ϕ]Lame[{b,-I b} Exp[-I ϕ]], r, Simplify]
T[{b,-I b} Exp[-I ϕ]]

κ0 = (λ+3μ)/ (λ+μ)

ClearAll[soln0];soln0[n_]:=Which[
n==0,{{r, 0},{0,r}}, 
n==1, {{1,I},r^2 {(κ0-2),-I (κ0+2)}}, 
n==-1,{{1,-I},r^2{(κ0-2),I (κ0+2)}}, 
n≥2, {n r^(n-1){1,I},r^(n+1) {(κ0-n -1),-I (κ0+n+1)}},
n≤-2, {Abs[n] r^(Abs[n]-1){1,-I },r^(Abs[n]+1) {(κ0-Abs[n]-1),I (κ0+Abs[n]+1)}}];

Formula for |n|>=2 in fact works for |n|=1, and we therefore redefine soln0:

({n r^(n-1){1,I},r^(n+1) {(κ0-n -1),-I (κ0+n+1)}}/.n->1)- {{1,I},r^2 {(κ0-2),-I (κ0+2)}}

soln0[n_]:=Piecewise[{
{{{r, 0},{0,r}}, n==0},
{{n r^(n-1){1,I},r^(n+1) {(κ0-n -1),-I (κ0+n+1)}}Exp[I n ϕ],n≥1},
{{Abs[n] r^(Abs[n]-1){1,-I },r^(Abs[n]+1) {(κ0-Abs[n]-1),I (κ0+Abs[n]+1)}}Exp[I n ϕ],n≤-1}
}];

Table[Map[Lame,soln0[n]], {n,-4,4}]//Simplify

Table[Map[T,soln0[n] ], {n,-2,2}]//Simplify

Table[Map[TDir,soln0[n]], {n,-2,2}]//Simplify

Solutions of the Lame system for Λ≠0. 
We seek these solutions in the form Grad[Φ[r, ϕ]] + Curl[(OverVector[z])Ψ[r,ϕ]]. Then the scalar potentials Φ and Ψ should satisfy Helmholtz equations, as seen from the calculation below. Indeed, setting ug=Grad[Φ[r, ϕ]] and uc=Curl[(OverVector[z])Ψ[r,ϕ]], we get
Lame[ug] - Λ ug = Grad[-(λ+2μ)ΔΦ - ΛΦ] 
and 
Lame[uc] - Λ uc = Curl[(OverVector[z])(-μΔΨ - ΛΨ)]

Curl[ Curl[{0, 0, Ψ[r, ϕ]}, {r, ϕ, z}, "Cylindrical"],{r, ϕ, z}, "Cylindrical"][[3]]
Laplacian[Ψ[r, ϕ], {r, ϕ, z}, "Cylindrical"]

ug = Grad[Φ[r, ϕ], {r, ϕ, z}, "Cylindrical"][[1 ;; 2]]
eqg1 = FullSimplify[Lame[ug] - Λ ug];
eqg2 = FullSimplify[
   Grad[-(λ+2μ)Laplacian[Φ[r, ϕ], {r, ϕ, z}, "Cylindrical"] - 
      Λ Φ[r, ϕ], {r, ϕ, z}, "Cylindrical"][[1 ;; 2]]];
eqg1 - eqg2

uc = Curl[{0, 0, Ψ[r, ϕ]}, {r, ϕ, z}, "Cylindrical"][[1 ;; 2]]
eqc1 = FullSimplify[Lame[uc] - Λ uc];
eqc2 = FullSimplify[
   Curl[{0, 0, -μ Laplacian[Ψ[r, ϕ], {r, ϕ, z}, "Cylindrical"] - Λ Ψ[r, ϕ]}, 
     {r, ϕ, z}, "Cylindrical"][[1 ;; 2]]];
eqc1 - eqc2//Simplify

The solutions of Helmholtz equation -ΔΘ-kΘ=0 are θ=(Subscript[J, |n|])((Sqrt[k])r)exp(I n ϕ) for k>0 and  θ=(Subscript[I, |n|])((Sqrt[-k])r)exp(I n ϕ)  for k<0. 
We will only deal with non-negative n’s remembering later to double the multiplicities for n≠0. solnneg refers to k<0 and may be ignored

solnpos[n_]:={Grad[BesselJ[n,Sqrt[Λ/(λ+2μ)]r] Exp[I n ϕ],{r,ϕ,z},"Cylindrical"][[1;;2]],Curl[{0,0,BesselJ[n,Sqrt[Λ/μ]r] Exp[I n ϕ]},{r,ϕ,z},"Cylindrical"][[1;;2]]}/.BesselJ[-1+n,z_]->(2 n BesselJ[n,z]/z-BesselJ[n+1,z])
solnneg[n_]:={Grad[BesselI[n,Sqrt[-Λ/(λ+2μ)]r] Exp[I n ϕ],{r,ϕ,z},"Cylindrical"][[1;;2]],Curl[{0,0,BesselI[n,Sqrt[-Λ/μ]r] Exp[I n ϕ]},{r,ϕ,z},"Cylindrical"][[1;;2]]}/.BesselI[-1+n,z_]->(2 n BesselI[n,z]/z+BesselI[n+1,z])//FullSimplify;

Just to check that we indeed produce solutions of Lame system for both positive and negative n:

Map[(Lame[#]-Λ #)&, solnpos[n]]//FullSimplify

solnposMinusn[n_] := {Grad[BesselJ[n,Sqrt[Λ/(λ+2μ)]r] Exp[-I n ϕ],{r,ϕ,z},"Cylindrical"][[1;;2]],Curl[{0,0,BesselJ[n,Sqrt[Λ/μ]r] Exp[-I n ϕ]},{r,ϕ,z},"Cylindrical"][[1;;2]]}/.BesselJ[-1+n,z_]->(2 n BesselJ[n,z]/z-BesselJ[n+1,z]);

Map[(Lame[#]-Λ #)&, solnposMinusn[n]]//FullSimplify

This is to check that choosing plus or minus n does not affect b.c. : the first (gradient) solutions changes the radial part as (u,v)->(u,-v) under n->-n,
the second (curl) solution changes the radial part as (u,v)->(-u,v) under n->-n, so homogeneous conditions stay the same. So from now on we just deal with n>= 0 and double all the eigenvalues for non-zero n

Exp[-I n ϕ] solnpos[n][[1]]-
Exp[I n ϕ]{{1,0},{0,-1}} .solnposMinusn[n][[1]]//Simplify
Exp[-I n ϕ] solnpos[n][[2]]-
Exp[I n ϕ]{{-1,0},{0,1}} .solnposMinusn[n][[2]]//Simplify
Exp[-I n ϕ] T[solnpos[n][[1]]]-
Exp[I n ϕ]{{1,0},{0,-1}} .T[solnposMinusn[n][[1]]]//Simplify
Exp[-I n ϕ] T[solnpos[n][[2]]]-
Exp[I n ϕ]{{-1,0},{0,1}} .T[solnposMinusn[n][[2]]]//Simplify

The Dirichlet eigenvalues are the Λ-roots of

(Map[(#/.{r->1, Exp[_]->1}) &, solnpos[n]]//Det)//FullSimplify

Direigseq=(Map[(#/.{r->1, Exp[_]->1}) &, solnpos[n]]//Det)/.BesselJ[-1+n,z_]->(2 n BesselJ[n,z]/z-BesselJ[n+1,z])//FullSimplify

Direigseqω=(Direigseq/.{(Λ/(λ+2 μ))->(Subscript[ω, 1,Λ]),(Λ/μ)->(Subscript[ω, 2,Λ]),n->k})
Direigseqω/.n->0
Direigseqω//TeXForm

and the Neumann eigenvalues are Λ-roots of (this sometimes does not show as Bessels but as Hypergeometric functions)

Neueigseq=(Map[(T[#]/.{r->1, Exp[_]->1}) &, solnpos[n]]//Det)/.BesselJ[-1+n,z_]->(2 n BesselJ[n,z]/z-BesselJ[n+1,z])//FullSimplify

Direq[λ_,μ_,n_,Λ_]:=n (Sqrt[(1/μ)]) BesselJ[n,(Sqrt[(Λ/(λ+2 μ))])] BesselJ[1+n,(Sqrt[(Λ/μ)])]+(Sqrt[(1/(λ+2 μ))]) (n BesselJ[n,(Sqrt[(Λ/μ)])]-(Sqrt[(Λ/μ)]) BesselJ[1+n,(Sqrt[(Λ/μ)])]) BesselJ[1+n,(Sqrt[(Λ/(λ+2 μ))])];

Simplify[Direq[λ,μ,n,Λ]-Direigseq/Sqrt[Λ], Λ>0&&μ>0&&λ+2 μ>0]

Direq[λ,μ,n,Λ]//FullSimplify
Direq[λ,μ,0,Λ]//Simplify

Neueq[λ_,μ_,n_,Λ_]:=2 (Sqrt[(Λ/μ)]) μ BesselJ[1+n,(Sqrt[(Λ/μ)])] ((Λ+2 n μ-2 (n^3) μ) BesselJ[n,(Sqrt[(Λ/(λ+2 μ))])]+2 (-1+(n^2)) μ (Sqrt[(Λ/(λ+2 μ))]) BesselJ[1+n,(Sqrt[(Λ/(λ+2 μ))])])+BesselJ[n,(Sqrt[(Λ/μ)])] (-Λ (Λ-4 (-1+n) n μ) BesselJ[n,(Sqrt[(Λ/(λ+2 μ))])]+2 μ (Sqrt[(Λ/(λ+2 μ))]) (Λ+2 n μ-2 (n^3) μ) BesselJ[1+n,(Sqrt[(Λ/(λ+2 μ))])]);

Simplify[Neueq[λ,μ,n,Λ]-Neueigseq, Λ>0&&μ>0&&λ+2 μ>0]

Direigs[λ0_,μ0_,nmax_, Λmax_]:=Module[{t,j},
t=Table[findAllRoots[Direq[λ0,μ0,n,Λ],{Λ,0.,Λmax}], {n,0,nmax}];
t={t[[1]], Table[{t[[j]],t[[j]]}, {j,2,Length[t]}]}//Flatten;
t//Sort
];
DireigsWithn[λ0_,μ0_, nmax_, Λmax_]:=Module[{t,t1,t2,t3,pair},
t=Table[findAllRoots[Direq[λ0,μ0,n,Λ],{Λ,0.,Λmax}], {n,0,nmax}];
t1=Select[Table[If[Length[t[[n+1]]]>0,Thread[{n,t[[n+1]]}]], {n,0,nmax}], Length[#]>0 &];
t2=SortBy[Flatten[t1,1],N[#[[2]]]&];
t2
];
Neueigs[λ0_,μ0_, nmax_, Λmax_]:=Module[{t,k},
t=Table[findAllRoots[Neueq[λ0,μ0,n,Λ],{Λ,0.,Λmax}], {n,0,nmax}];
t={t[[1]],Table[{t[[j]],t[[j]]}, {j,2,Length[t]}]}//Flatten;
{0.,0.,0.,t//Sort}//Flatten 
]
NeueigsWithn[λ0_,μ0_, nmax_, Λmax_]:=Module[{t,t1,t2,t3,pair},
t=Table[findAllRoots[Neueq[λ0,μ0,n,Λ],{Λ,0.,Λmax}], {n,0,nmax}];
t1=Select[Table[If[Length[t[[n+1]]]>0,Thread[{n,t[[n+1]]}]], {n,0,nmax}], Length[#]>0 &];
t2=SortBy[Flatten[t1,1],#[[2]]&];
PrependTo[t2,{1,0.}];
PrependTo[t2,{0,0.}];
t2
];

Examples:

Direigs[0,1, 52, 3200.]
DireigsWithn[0,1, 52, 3200.]

Neueigs[0,1, 52, 3200.]
NeueigsWithn[0,1, 52, 3200.]

MyNlambdaDisk[evs_,Λ_]:=Sum[If[evs[[j]][[1]]==0,1,2]UnitStep[Λ-evs[[j]][[2]]],{j,1,Length[evs]}];

Finding the number of harmonics to stop at

findAllRoots[Direq[-0.5,1,47,Λ],{Λ,0.,3000}]
findAllRoots[Direq[-0.5,1,48,Λ],{Λ,0.,3000}]

findAllRoots[Direq[3,1,46,Λ],{Λ,0.,3000}]
findAllRoots[Direq[3,1,47,Λ],{Λ,0.,3000}]

findAllRoots[Direq[100,1,46,Λ],{Λ,0.,3000}]
findAllRoots[Direq[100,1,47,Λ],{Λ,0.,3000}]

Dirichlet counting functions

evD1=DireigsWithn[-0.5,1, 47, 3000];
evD2=DireigsWithn[3,1, 46, 3000];
evD3=DireigsWithn[100,1, 46, 3000];

xtdisk=Table[{x, MaTeX[x]}, {x,1000,3000, 1000}];
ytdisk=Table[{x, MaTeX[x]}, {x,-60,-20, 20}];

figdisk1=Plot[{MyNlambdaDisk[evD1,Λ]-CWeyl[2,-0.5, 1]Pi Λ,2Pi BDir[2, αλ[-0.5,1],1] Sqrt[Λ], 2Pi BDirLiu[2,-0.5,1]Sqrt[Λ]},{Λ,0,3000},PlotStyle->{Black, Blue, Red},Exclusions->None,Ticks->{xtdisk,ytdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=-0.5"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,0}],Scaled[{0,0}]]];
figdisk2=Plot[{MyNlambdaDisk[evD2,Λ]-CWeyl[2,3, 1]Pi Λ,2Pi BDir[2, αλ[3,1],1] Sqrt[Λ], 2Pi BDirLiu[2,3,1]Sqrt[Λ]},{Λ,0,3000},PlotStyle->{Black, Blue, Red},Exclusions->None, Ticks->{xtdisk,ytdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=3"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,0}],Scaled[{0,0}]]]; 
figdisk3=Plot[{MyNlambdaDisk[evD3,Λ]-CWeyl[2,100, 1]Pi Λ,2Pi BDir[2, αλ[100,1],1] Sqrt[Λ], 2Pi BDirLiu[2,100,1]Sqrt[Λ]},{Λ,0,3000},PlotStyle->{Black, Blue, Red},Exclusions->None,Ticks->{xtdisk,ytdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=100"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,0}],Scaled[{0,0}]]] ;
legendisk=LineLegend[{Black, Blue, Red},MaTeX[{"\\text{computed }\\mathcal{N}_{\\mathrm{Dir}}(\\Lambda)-a\\pi \\Lambda\\qquad\\qquad","2\\pi b_{\\mathrm{Dir}}\\sqrt{\\Lambda}\\qquad\\qquad","2\\pi b_{\\mathrm{Dir}}^{\\mathrm{Liu}}\\sqrt{\\Lambda}"}],LegendFunction->None,LegendLayout->{"Row",1}];
figdisk=GraphicsColumn[{GraphicsRow[{figdisk1,figdisk2,figdisk3},ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
legendisk},ImageSize->Full,Spacings->{0,Automatic}]

Export["figdisk.pdf",figdisk]

Finding the number of harmonics to stop at

findAllRoots[Neueq[-0.5,1,56,Λ],{Λ,0.,2000}]
findAllRoots[Neueq[-0.5,1,57,Λ],{Λ,0.,2000}]

findAllRoots[Neueq[3,1,45,Λ],{Λ,0.,2000}]
findAllRoots[Neueq[3,1,46,Λ],{Λ,0.,2000}]

findAllRoots[Neueq[100,1,44,Λ],{Λ,0.,2000}]
findAllRoots[Neueq[100,1,45,Λ],{Λ,0.,2000}]

evDfr1=NeueigsWithn[-0.5,1, 56, 2000];
evDfr2=NeueigsWithn[3,1, 45, 2000];
evDfr3=NeueigsWithn[100,1, 44, 2000];

ytfrdisk=Table[{x, MaTeX["\\ "<>ToString[x]]}, {x,20,60, 20}];

figdiskfr1=Plot[{MyNlambdaDisk[evDfr1,Λ]-CWeyl[2,-0.5, 1]Pi Λ,2Pi Bfree[2, αλ[-0.5,1],1] Sqrt[Λ]},{Λ,0,2000},PlotStyle->{Black, Blue},Exclusions->None,Ticks->{xtdisk,ytfrdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=-0.5"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]];
figdiskfr2=Plot[{MyNlambdaDisk[evDfr2,Λ]-CWeyl[2,3, 1]Pi Λ,2Pi Bfree[2, αλ[3,1],1] Sqrt[Λ]},{Λ,0,2000},PlotStyle->{Black, Blue},Exclusions->None,Ticks->{xtdisk,ytfrdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=3"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]]; 
figdiskfr3=Plot[{MyNlambdaDisk[evDfr3,Λ]-CWeyl[2,100, 1]Pi Λ,2Pi Bfree[2, αλ[100,1],1] Sqrt[Λ]},{Λ,0,2000},PlotStyle->{Black, Blue},Exclusions->None,Ticks->{xtdisk,ytfrdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=100"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]] ;
legendiskfr=LineLegend[{Black, Blue},MaTeX[{"\\text{computed }\\mathcal{N}_{\\mathrm{free}}(\\Lambda)-a\\pi \\Lambda\\qquad\\qquad","2\\pi b_{\\mathrm{free}}\\sqrt{\\Lambda}\\qquad\\qquad"}],LegendFunction->None,LegendLayout->{"Row",1}];
figdiskfr=GraphicsColumn[{GraphicsRow[{figdiskfr1,figdiskfr2,figdiskfr3},ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
legendiskfr},ImageSize->Full,Spacings->{0,Automatic}]

Export["figdiskfree.pdf",figdiskfr]

figdiskfr1new=Plot[{MyNlambdaDisk[evDfr1,Λ]-CWeyl[2,-0.5, 1]Pi Λ,2Pi Bfree[2, αλ[-0.5,1],1] Sqrt[Λ],-2Pi BDirLiu[2, -0.5,1] Sqrt[Λ]},{Λ,0,2000},PlotStyle->{Black, Blue, Red},Exclusions->None,Ticks->{xtdisk,ytfrdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=-0.5"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]];
figdiskfr2new=Plot[{MyNlambdaDisk[evDfr2,Λ]-CWeyl[2,3, 1]Pi Λ,2Pi Bfree[2, αλ[3,1],1] Sqrt[Λ],-2Pi BDirLiu[2, 3,1] Sqrt[Λ]},{Λ,0,2000},PlotStyle->{Black, Blue,Red},Exclusions->None,Ticks->{xtdisk,ytfrdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=3"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]]; 
figdiskfr3new=Plot[{MyNlambdaDisk[evDfr3,Λ]-CWeyl[2,100, 1]Pi Λ,2Pi Bfree[2, αλ[100,1],1] Sqrt[Λ],-2Pi BDirLiu[2, 100,1] Sqrt[Λ]},{Λ,0,2000},PlotStyle->{Black, Blue,Red},Exclusions->None,Ticks->{xtdisk,ytfrdisk}, AxesLabel->{MaTeX["\\Lambda"],None},Epilog->Inset[Framed[MaTeX["\\lambda=100"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]] ;
legendiskfrnew=LineLegend[{Black, Blue,Red},MaTeX[{"\\text{computed }\\mathcal{N}_{\\mathrm{free}}(\\Lambda)-a\\pi \\Lambda\\qquad\\qquad","2\\pi b_{\\mathrm{free}}\\sqrt{\\Lambda}\\qquad\\qquad","2\\pi b_{\\mathrm{free}}^{\\mathrm{Liu}}\\sqrt{\\Lambda}"}],LegendFunction->None,LegendLayout->{"Row",1}];
figdiskfrnew=GraphicsColumn[{GraphicsRow[{figdiskfr1new,figdiskfr2new,figdiskfr3new},ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
legendiskfrnew},ImageSize->Full,Spacings->{0,Automatic}]

Export["figdiskfreenew.pdf",figdiskfrnew]