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

αλ[λ_,μ_]:=μ/(λ+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μ)μ(λ+μ));

Clear[findAllRoots]
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->{}]

ClearAll[ℒ];ℒ[V_, λ_, μ_]:= μ Curl[Curl[V,{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}]-(λ+2μ) Grad[Div[V,{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}];

ℒ[{(Subscript[v, 1])[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])],(Subscript[v, 2])[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])],(Subscript[v, 3])[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]}, λ, μ]

1/(λ+2 μ) ℒ[Grad[ψ[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}], λ, μ]+Laplacian[Grad[ψ[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}] 
1/μ ℒ[Curl[{0,0,ψ[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}], λ, μ] + Laplacian[Curl[{0,0,ψ[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}], {(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}]
1/μ ℒ[Curl[Curl[{0,0,ψ[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}], λ, μ] + Laplacian[Curl[Curl[{0,0,ψ[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}], {(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}]

ClearAll[κ];-Laplacian[Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])-Sqrt[κ-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))], {(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}]/Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])-Sqrt[κ-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))]//Simplify

V=Grad[(Subscript[c, 1])Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])+Sqrt[Λ/(λ+2 μ)-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))]   +(Subscript[c, 2])Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])-Sqrt[Λ/(λ+2 μ)-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}];
W=Curl[{0,0,(Subscript[c, 3])Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])+Sqrt[Λ/μ-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))] +(Subscript[c, 4])Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])-Sqrt[Λ/ μ-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))]},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}];
Y=Curl[Curl[{0,0,(Subscript[c, 5])Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])+Sqrt[Λ/μ-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))] +(Subscript[c, 6])Exp[I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2])-Sqrt[Λ/ μ-(Subscript[ξ, 1])^2-(Subscript[ξ, 2])^2] (Subscript[x, 3]))]},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}],{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}];
ℒ[V,  λ, μ]-Λ V //Simplify
ℒ[W,  λ, μ]-Λ W //Simplify
ℒ[Y,  λ, μ]-Λ Y //Simplify

U=V+W+Y

U0=(U/.{(Subscript[ξ, 1])->0, (Subscript[ξ, 2])->0})
CoefficientArrays[{U0[[3]]/.(Subscript[x, 3])->0, U0[[3]]/.(Subscript[x, 3])->h}, {(Subscript[c, 1]), (Subscript[c, 2])}][[2]]//Normal//Det//FullSimplify

f=((CoefficientArrays[{(Exp[-I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2]))] U/.(Subscript[x, 3])->0),(Exp[-I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2]))] U/.(Subscript[x, 3])->h)}//Flatten//Simplify, Table[(Subscript[c, j]), {j,1,6}]][[2]]//Normal//Det//Simplify)/.{(Subscript[ξ, 1])->m1,(Subscript[ξ, 2])->m2, m1^2+m2^2->R,- m1^2-m2^2->-R} )//FullSimplify;
f1=f/.{m1^2+m2^2->R,- m1^2-m2^2->-R}//Simplify;
mydet1=f1/.{(Sqrt[-R+(Λ/(λ+2 μ))])->Sqrt[C1],(( (Sqrt[Λ-R μ]))/(Sqrt[μ]))->Sqrt[C2], (Sqrt[Λ-R μ])-> (Sqrt[C2])Sqrt[μ]}//FullSimplify

(* Case Λ< R μ    *)
(mydet1/.{Sqrt[C1]-> I Sqrt[Abs[C1]], Sqrt[C2]-> I Sqrt[Abs[C2]]})//Simplify

(* Case   R μ <= Λ < R (λ+2 μ) *)
(mydet1/.{Sqrt[C1]-> I Sqrt[Abs[C1]]})//Simplify

MyDetSign1[λ_, μ_, h_, R_ ,Λ_]:=  (* for Λ < R μ *) With[{AC1=R-(Λ/(λ+2 μ)),AC2= -Λ/μ + R}, (2 R (Sqrt[AC1]) (Sqrt[AC2]) (-1+Cosh[h (Sqrt[AC1])] Cosh[h (Sqrt[AC2])])-((((Λ^2)+2 (R^2) μ (λ+2 μ)-R Λ (λ+3 μ)) Sinh[h (Sqrt[AC1])] Sinh[h (Sqrt[AC2])])/(μ (λ+2 μ))))/(Cosh[h (Sqrt[AC1])] Cosh[h (Sqrt[AC2])])];
MyDetSign2Part1[λ_, μ_, h_, R_ ,Λ_]:=  (* for R μ < Λ < R(λ+2 μ)  *) With[{AC1=R-(Λ/(λ+2 μ)),AC2= Λ/ μ-R}, Sin[(Sqrt[AC2]) h]];
MyDetSign2Part2[λ_, μ_, h_, R_ ,Λ_]:=  (* for R μ < Λ < R(λ+2 μ)  *) With[{AC1=R-(Λ/(λ+2 μ)),AC2= Λ/ μ-R}, (-2 (-1+Cosh[h (Sqrt[AC1])] Cos[h (Sqrt[AC2])])+((((Λ^2)+2 (R^2) μ (λ+2 μ)-R Λ (λ+3 μ)) Sinh[h (Sqrt[AC1])] Sin[h (Sqrt[AC2])])/(μ (λ+2 μ)R (Sqrt[AC1]) (Sqrt[AC2]))))/Cosh[h (Sqrt[AC1])] ];
MyDetSign3Part1[λ_, μ_, h_, R_ ,Λ_]:=  (* for Λ > R(λ+2 μ)  *) With[{C1=-R+(Λ/(λ+2 μ)),C2= Λ/ μ-R}, 8Sin[(Sqrt[C2]) h]];
MyDetSign3Part2[λ_, μ_, h_, R_ ,Λ_]:=  (* for Λ > R(λ+2 μ)  *) With[{C1=-R+(Λ/(λ+2 μ)),C2= Λ/ μ-R},
(2  (-1+Cos[(Sqrt[C1]) h] Cos[(Sqrt[C2]) h])-((((Λ^2)+2 (R^2) μ (λ+2 μ)-R Λ (λ+3 μ)) Sin[(Sqrt[C1]) h] Sin[(Sqrt[C2]) h])/(μ (λ+2 μ)(Sqrt[C1]) (Sqrt[C2]) R))) ];

(*Example:λ=0,μ=1,h=4,R=9 *)
Plot[MyDetSign1[0, 1, 4, 9,Λ], {Λ, 0, 9}]
Show[Plot[{MyDetSign2Part1[0, 1, 4, 9,Λ],MyDetSign2Part2[0, 1, 4, 9,Λ]}, {Λ, 9, 18}],
ListPlot[Thread[{Λ/.NSolve[MyDetSign2Part2[0, 1, 4, 9,Λ]==0 && 9<Λ<18],0}],PlotStyle->{{Orange,PointSize[Medium]}}]]
Show[Plot[{MyDetSign3Part1[0, 1, 4, 9,Λ],MyDetSign3Part2[0, 1, 4, 9,Λ]}, {Λ, 18, 100}],
ListPlot[Thread[{Λ/.NSolve[MyDetSign3Part2[0, 1, 4, 9,Λ]==0 && 18<Λ<100],0}],PlotStyle->{{Orange,PointSize[Medium]}}]]

(* numerically computed eigenvalues *)
MyEigs[λ_, μ_, h_,ΛMax_]:=Module[{evs, R, Rmax, n, Rsq, solns21,solns22, solns31,solns32, j},
evs=Table[{0, 1,N[(λ+2 μ)Pi^2 n^2/h^2]},{n,1,h Sqrt[ΛMax]/(Pi Sqrt[λ+2 μ])}];
(* the eigenvalues above come from R=0, with multiplicity one *)
Do[Rsq=SquaresR[2,R];
If[Rsq!=0,
solns21=Λ/.NSolve[MyDetSign2Part1[λ,μ,h,R,Λ]==0 && μ R<Λ<(λ+2μ)R && Λ<=ΛMax, Λ];
solns22=Λ/.NSolve[MyDetSign2Part2[λ,μ,h,R,Λ]==0 && μ R<Λ<(λ+2μ)R && Λ<=ΛMax, Λ];
Do[AppendTo[evs,{R,Rsq, solns21[[j]]}],{j,1,Length[solns21]}];
Do[AppendTo[evs,{R,Rsq, solns22[[j]]}],{j,1,Length[solns22]}];
solns31=Λ/.NSolve[MyDetSign3Part1[λ,μ,h,R,Λ]==0 &&(λ+2μ)R < Λ<=ΛMax, Λ];
solns32=Λ/.NSolve[MyDetSign3Part2[λ,μ,h,R,Λ]==0 &&(λ+2μ)R < Λ<=ΛMax, Λ];
Do[AppendTo[evs,{R,Rsq, solns31[[j]]}],{j,1,Length[solns31]}];
Do[AppendTo[evs,{R,Rsq, solns32[[j]]}],{j,1,Length[solns32]}];
]
,{R,1,ΛMax/μ}];
SortBy[evs, N[#[[3]]] &]
];
MyEigs3D2D[evs_]:=Module[{evs2,j,Rsq2},
evs2={};
Do[Rsq2=SquaresR[1,evs[[j]][[1]]];If[Rsq2!=0,AppendTo[evs2,{evs[[j]][[1]],Rsq2, evs[[j]][[3]]}]],{j,1,Length[evs]}];
evs2
]

λs={-1/2,3,100}; hs={1/2,2}; Lmax={1000,500,250};
htext={"0.5", "2"}; ylim={-5000, -2500, -1200}; lambdatext={"-0.5", "3", "100"};
xticks=Table[{{Lmax[[j]]/2, MaTeX[Lmax[[j]]/2]},{Lmax[[j]], MaTeX[Lmax[[j]]]}},{j,1,3}] ;
yticks=Table[{{ylim[[i]]/2, MaTeX[ylim[[i]]/2]},{ylim[[i]], MaTeX[ylim[[i]]]}},{i,1,2}] ;

EVs=Table[MyEigs[λs[[j]],1,hs[[i]],Lmax[[i]]],{i,1,2},{j,1,3}];

MyNlambda3d[evs_,Λ_]:=Sum[evs[[j]][[2]]UnitStep[Λ-evs[[j]][[3]]],{j,1,Length[evs]}];
Vol3[h_]:=(2Pi)^2 h;
Area3[h_]:=2(2Pi)^2;

graphicstableDir=
Table[
Plot[{MyNlambda3d[EVs[[i]][[j]],Λ]-CWeyl[3,λs[[j]],1]Vol3[hs[[i]]]  Λ^(3/2), BDir[3,αλ[λs[[j]],1],1]Area3[hs[[i]]]Λ, BDirLiu[3,λs[[j]],1]Area3[hs[[i]]]Λ},{Λ, 0, Lmax[[i]]}, PlotStyle->{Black,Blue,Red},
Axes->True,
AxesLabel->{MaTeX["\\Lambda"], None},
PlotRange->{{0,Lmax[[i]]},{ylim[[i]],0}},
Ticks->{xticks[[i]],yticks[[i]]},
Epilog->Inset[Framed[MaTeX["\\lambda="<>lambdatext[[j]]<>",\\quad h="<>htext[[i]]], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,0}],Scaled[{0,0}]]
],{i,1,2},{j,1,3}]
(* plotsCylinder3=GraphicsGrid[graphicstableDir, ImageSize->Full, Frame->True, Dividers->All] *)

legencyl=LineLegend[{Black, Blue, Red},MaTeX[{"\\text{computed }\\mathcal{N}_{\\mathrm{Dir}}(\\Lambda)- a (2\\pi)^2 h \\Lambda^{3/2}\\qquad\\qquad","8\\pi^2 b_{\\mathrm{Dir}}\\Lambda\\qquad\\qquad","8\\pi^2 b_{\\mathrm{Dir}}^{\\mathrm{Liu}}\\Lambda"}],LegendFunction->None,LegendLayout->{"Row",1}]

figcyl=GraphicsColumn[{GraphicsRow[graphicstableDir[[1]],ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
GraphicsRow[graphicstableDir[[2]],ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
legencyl},ImageSize->Full,Spacings->{0,Automatic}]

Export["figcyl.pdf",figcyl]

TU= {μ (D[U[[1]], (Subscript[x, 3])]+D[U[[3]], (Subscript[x, 1])]),μ(D[U[[2]], (Subscript[x, 3])]+D[U[[3]], (Subscript[x, 2])]),λ(D[U[[1]], (Subscript[x, 1])]+D[U[[2]], (Subscript[x, 2])]+D[U[[3]], (Subscript[x, 3])])+2μ D[U[[3]], (Subscript[x, 3])]}Exp[-I ((Subscript[ξ, 1]) (Subscript[x, 1])+(Subscript[ξ, 2]) (Subscript[x, 2]))]//Simplify;

TU0={μ (D[U0[[1]], (Subscript[x, 3])]+D[U0[[3]], (Subscript[x, 1])]),μ(D[U0[[2]], (Subscript[x, 3])]+D[U0[[3]], (Subscript[x, 2])]),λ(D[U0[[1]], (Subscript[x, 1])]+D[U0[[2]], (Subscript[x, 2])]+D[U0[[3]], (Subscript[x, 3])])+2μ D[U0[[3]], (Subscript[x, 3])]}//Simplify;
CoefficientArrays[{TU0[[3]]/.(Subscript[x, 3])->0, TU0[[3]]/.(Subscript[x, 3])->h}, {(Subscript[c, 1]), (Subscript[c, 2])}][[2]]//Normal//Det//FullSimplify

Tfmatrix=(((CoefficientArrays[{(TU/.(Subscript[x, 3])->0),(TU/.(Subscript[x, 3])->h)}//Flatten//Simplify, Table[(Subscript[c, j]), {j,1,6}]][[2]]//Normal)/.{(Subscript[ξ, 1])->m1,(Subscript[ξ, 2])->m2, m2^2->R-m1^2,-m2^2->-R+m1^2} )//Simplify)/.{(Subscript[ξ, 1])->m1,(Subscript[ξ, 2])->m2, m2^2->R-m1^2,-m2^2->-R+m1^2} //Simplify

Tmatrix1=(((Tfmatrix/.Λ->μ R)//FullSimplify)/.(Sqrt[-((R (λ+μ))/(λ+2 μ))])->I(Sqrt[((R (λ+μ))/(λ+2 μ))]))//Simplify;
Tmatrix1//TableForm
Tmatrix1//NullSpace

U/.{(Subscript[c, 1])->0, (Subscript[c, 2])->0,(Subscript[c, 3])->0,(Subscript[c, 4])->0, (Subscript[c, 5])->-1,(Subscript[c, 6])->1, (Λ/μ)-(Subscript[ξ, 1]^2)-(Subscript[ξ, 2]^2)->0}//Simplify
U1=U/.{(Subscript[c, 1])->0, (Subscript[c, 2])->0,(Subscript[c, 3])->1,(Subscript[c, 4])->0, (Subscript[c, 5])->0,(Subscript[c, 6])->0, (Λ/μ)-(Subscript[ξ, 1]^2)-(Subscript[ξ, 2]^2)->0}//FullSimplify
U/.{(Subscript[c, 1])->0, (Subscript[c, 2])->0,(Subscript[c, 3])->0,(Subscript[c, 4])->1, (Subscript[c, 5])->0,(Subscript[c, 6])->0, (Λ/μ)-(Subscript[ξ, 1]^2)-(Subscript[ξ, 2]^2)->0}//FullSimplify

ℒ[U1,  λ, μ][[1;;2]]/U1[[1;;2]]//Simplify

U11=ComplexExpand[U1+Conjugate[U1]]/2//Factor
U12=ComplexExpand[U1-Conjugate[U1]]/2//Factor

{U1/.{(Subscript[ξ, 1])->2, (Subscript[ξ, 2])->0},U1/.{(Subscript[ξ, 1])->-2, (Subscript[ξ, 2])->0},U1/.{(Subscript[ξ, 2])->2, (Subscript[ξ, 1])->0},U1/.{(Subscript[ξ, 2])->-2, (Subscript[ξ, 1])->0}}//Simplify

{U1/.{(Subscript[ξ, 1])->2, (Subscript[ξ, 2])->1},U1/.{(Subscript[ξ, 1])->-2, (Subscript[ξ, 2])->-1},U1/.{(Subscript[ξ, 1])->2, (Subscript[ξ, 2])->-1},U1/.{(Subscript[ξ, 1])->-2, (Subscript[ξ, 2])->1},U1/.{(Subscript[ξ, 2])->2, (Subscript[ξ, 1])->1},U1/.{(Subscript[ξ, 2])->-2, (Subscript[ξ, 1])->-1},U1/.{(Subscript[ξ, 2])->2, (Subscript[ξ, 1])->-1},U1/.{(Subscript[ξ, 2])->-2, (Subscript[ξ, 1])->1}}//Simplify//ComplexExpand

Tmatrix2=(((Tfmatrix/.Λ->(λ+2μ) R)//FullSimplify))//Simplify;
Tmatrix2//TableForm
Tmatrix2//NullSpace

U/.{(Subscript[c, 1])->-1, (Subscript[c, 2])->1,(Subscript[c, 3])->0,(Subscript[c, 4])->0, (Subscript[c, 5])->0,(Subscript[c, 6])->0,Λ ->(λ+2μ)((Subscript[ξ, 1]^2)+(Subscript[ξ, 2]^2))}//FullSimplify

Tf=((CoefficientArrays[{(TU/.(Subscript[x, 3])->0),(TU/.(Subscript[x, 3])->h)}//Flatten//Simplify, Table[(Subscript[c, j]), {j,1,6}]][[2]]//Normal//Det)/.{(Subscript[ξ, 1])->m1,(Subscript[ξ, 2])->m2, m2^2->R-m1^2,-m2^2->-R+m1^2} )//Simplify;

Tf1=(Tf/.{m2^2->R-m1^2,-m2^2->-R+m1^2,m2^4->R^2-2R m1^2+m1^4})//Simplify

mydetT1=(Tf1/.{(Sqrt[-R+(Λ/(λ+2 μ))])->Sqrt[C1],(( (Sqrt[Λ-R μ]))/(Sqrt[μ]))->Sqrt[C2], (Sqrt[-R+(Λ/μ)])->Sqrt[C2],(Sqrt[Λ-R μ])-> (Sqrt[C2])Sqrt[μ]})//Simplify

mydetT2=(mydetT1//ExpToTrig//TrigExpand//Factor)

mydetT3=mydetT2/((1/(λ+2 μ))4 I (Sqrt[C2]) (R^2) (μ^2) )//FullSimplify

mydetT3/.{C2->-AC2, C1->-AC1}/.{Sqrt[-AC2]->I Sqrt[AC2],  Sqrt[-AC1]->I Sqrt[AC1]}//Simplify
mydetT3/.{ C1->-AC1, C2->AC2}/.{Sqrt[-AC1]->I Sqrt[AC1]}//Factor//Simplify

MyDetTSign1[λ_, μ_, h_, R_ ,Λ_]:=  (* for Λ < R μ *) With[{AC1=R-(Λ/(λ+2 μ)),AC2= -Λ/μ + R}, (-8 (Sqrt[AC1]) R μ (λ+2 μ) (((Λ-2 R μ))^2) (Λ-R μ) (-1+Cosh[(Sqrt[AC1]) h] Cosh[(Sqrt[AC2]) h])-(Sqrt[AC2]) (λ (Λ^4)+2 (Λ^3) (-4 R λ+Λ) μ+8 R (3 R λ-2 Λ) (Λ^2) (μ^2)+16 (R^2) Λ (-3 R λ+4 Λ) (μ^3)+16 (R^3) (2 R λ-7 Λ) (μ^4)+64 (R^4) (μ^5)) Sinh[(Sqrt[AC1]) h] Sinh[(Sqrt[AC2]) h])/(Cosh[h (Sqrt[AC1])]Cosh[h (Sqrt[AC2])])];
MyDetTSign2Part1[λ_, μ_, h_, R_ ,Λ_]:=  (* for R μ < Λ < R(λ+2 μ)  *) With[{AC1=R-(Λ/(λ+2 μ)),AC2= Λ/ μ-R}, Sin[(Sqrt[AC2]) h]];
MyDetTSign2Part2[λ_, μ_, h_, R_ ,Λ_]:=  (* for R μ < Λ < R(λ+2 μ)  *) With[{AC1=R-(Λ/(λ+2 μ)),AC2= Λ/ μ-R}, (8 (Sqrt[AC1]) R μ (λ+2 μ) (((Λ-2 R μ))^2) (Λ-R μ)+8 (Sqrt[AC1]) R μ (λ+2 μ) (((Λ-2 R μ))^2) (-Λ+R μ) Cos[(Sqrt[AC2]) h] Cosh[(Sqrt[AC1]) h]+(Sqrt[AC2]) (2 μ ((Λ^4)-8 R (Λ^3) μ+32 (R^2) (Λ^2) (μ^2)-56 (R^3) Λ (μ^3)+32 (R^4) (μ^4))+λ ((Λ^4)-8 R (Λ^3) μ+24 (R^2) (Λ^2) (μ^2)-48 (R^3) Λ (μ^3)+32 (R^4) (μ^4))) Sin[(Sqrt[AC2]) h] Sinh[(Sqrt[AC1]) h])/( Cosh[h (Sqrt[AC1])]Λ)^2];
MyDetTSign3Part1[λ_, μ_, h_, R_ ,Λ_]:=  (* for Λ > R(λ+2 μ)  *) With[{C1=-R+(Λ/(λ+2 μ)),C2= Λ/ μ-R}, 2Sin[(Sqrt[C2]) h]];
MyDetTSign3Part2[λ_, μ_, h_, R_ ,Λ_]:=  (* for Λ > R(λ+2 μ)  *) With[{C1=-R+(Λ/(λ+2 μ)),C2= Λ/ μ-R},
(8 (Sqrt[C1]) R μ (λ+2 μ) (((Λ-2 R μ))^2) (Λ-R μ) (-1+Cos[(Sqrt[C1]) h] Cos[(Sqrt[C2]) h])-(Sqrt[C2]) (λ (Λ^4)+2 (Λ^3) (-4 R λ+Λ) μ+8 R (3 R λ-2 Λ) (Λ^2) (μ^2)+16 (R^2) Λ (-3 R λ+4 Λ) (μ^3)+16 (R^3) (2 R λ-7 Λ) (μ^4)+64 (R^4) (μ^5)) Sin[(Sqrt[C1]) h] Sin[(Sqrt[C2]) h]) ]/Λ^4;

(* numerically computed eigenvalues *)
MyTEigs[λ_, μ_, h_,ΛMax_]:=Module[{evs, R, Rmax, n, Rsq, solns1,solns21,solns22, solns31,solns32, j, k,  i6, M, M2},
PrintTemporary["λ=",λ," μ=",μ];
evs=Table[{0, 1,N[(λ+2 μ)Pi^2 n^2/h^2]},{n,1,h Sqrt[ΛMax]/(Pi Sqrt[λ+2 μ])}];
(* the eigenvalues above come from R=0, with multiplicity one *)
evs=AppendTo[evs,{0,4,0}]; (*ev=0*)
Monitor[Do[Rsq=SquaresR[2,R];
If[Rsq!=0,
M=Min[μ R, ΛMax]-10^(-6);
solns1=findAllRoots[MyDetTSign1[λ,μ,h,R,Λ]/R^3, {Λ,10^(-4),M}, PlotPoints->10000];
Do[AppendTo[evs,{R,Rsq, solns1[[j]]}],{j,1,Length[solns1]}];
AppendTo[evs, {R, Rsq, μ R}];
If[R<=ΛMax/μ,
M2=Min[(λ+2μ)R, ΛMax]-10^(-6);
solns21=findAllRoots[MyDetTSign2Part1[λ,μ,h,R,Λ],{Λ, μ R+10^(-4),M2}, PlotPoints->10000];
solns22=findAllRoots[MyDetTSign2Part2[λ,μ,h,R,Λ],{Λ, μ R+10^(-4),M2}, PlotPoints->10000];
Do[AppendTo[evs,{R,Rsq, solns21[[j]]}],{j,1,Length[solns21]}];
Do[AppendTo[evs,{R,Rsq, solns22[[j]]}],{j,1,Length[solns22]}];
solns31=findAllRoots[MyDetTSign3Part1[λ,μ,h,R,Λ],{Λ,(λ+2μ)R+10^(-4),ΛMax}, PlotPoints->10000];
solns32=findAllRoots[MyDetTSign3Part2[λ,μ,h,R,Λ],{Λ,(λ+2μ)R+10^(-4),ΛMax}, PlotPoints->10000];
Do[AppendTo[evs,{R,Rsq, solns31[[j]]}],{j,1,Length[solns31]}];
Do[AppendTo[evs,{R,Rsq, solns32[[j]]}],{j,1,Length[solns32]}];
];
]
,{R,1,3/2 ΛMax/μ}],ProgressIndicator[R, {1,3/2 ΛMax/μ}]];
SortBy[evs, N[#[[3]]] &]
];

λs={-1/2,3,100}; hs={1/2,2}; Lmax={1000,500,250};
htext={"0.5", "2"}; ylim={5000, 2500, 1200}; lambdatext={"-0.5", "3", "100"};
xticks=Table[{{Lmax[[j]]/2, MaTeX[Lmax[[j]]/2]},{Lmax[[j]], MaTeX[Lmax[[j]]]}},{j,1,3}] ;
yticks=Table[{{ylim[[i]]/2, MaTeX[ylim[[i]]/2]},{ylim[[i]], MaTeX[ylim[[i]]]}},{i,1,2}] ;

(* EVTs=Table[MyTEigs[λs[[j]],1,hs[[i]],Lmax[[i]]],{i,1,2},{j,1,3}]; *)

Lmax1={300,400,500};
EVTs1=Table[MyTEigs[λs[[j]],1,1,1.1Lmax1[[j]]],{j,1,3}];

(* graphicstableFree=
Table[
Plot[{MyNlambda3d[EVTs[[i]][[j]],Λ]-CWeyl[3,λs[[j]],1]Vol3[hs[[i]]]  Λ^(3/2), Bfree[3,αλ[λs[[j]],1],1]Area3[hs[[i]]]Λ},{Λ, 0, Lmax[[i]]}, PlotStyle->{Black,Blue}, Exclusions->None,
Axes->True,
AxesLabel->{MaTeX["\\Lambda"], None},
PlotRange->{{0,Lmax[[i]]},{0,ylim[[i]]}},
Ticks->{xticks[[i]],yticks[[i]]},
Epilog->Inset[Framed[MaTeX["\\lambda="<>lambdatext[[j]]<>", h="<>htext[[i]]], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]] 
],{i,1,2},{j,1,3}]  *)

xticks1=Table[{{Lmax1[[j]]/2, MaTeX[Lmax1[[j]]/2]},{Lmax1[[j]], MaTeX[Lmax1[[j]]]}},{j,1,3}] ;
yticks1={{1000, MaTeX[1000]}, {2000, MaTeX[2000]}};

graphicstable1Free=
Table[
Plot[{MyNlambda3d[EVTs1[[j]],Λ]-CWeyl[3,λs[[j]],1]Vol3[1]  Λ^(3/2), Bfree[3,αλ[λs[[j]],1],1]Area3[1]Λ},{Λ, 0, Lmax1[[j]]}, PlotStyle->{Black,Blue},Exclusions->None,
Axes->True,
AxesLabel->{MaTeX["\\Lambda"], None},
PlotRange->{{0,Lmax1[[j]]},Automatic},
Ticks->{xticks1[[j]],yticks1},
Epilog->Inset[Framed[MaTeX["\\lambda="<>lambdatext[[j]]<>",\\quad h=1"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]
],{j,1,3}]

legendcylfree=LineLegend[{Black, Blue},MaTeX[{"\\text{computed }\\mathcal{N}_{\\mathrm{free}}(\\Lambda)- a (2\\pi)^2 h \\Lambda^{3/2}\\qquad\\qquad","8\\pi^2 b_{\\mathrm{free}}\\Lambda"}],LegendFunction->None,LegendLayout->{"Row",1}]

graphicstable1Freenew=
Table[
Plot[{MyNlambda3d[EVTs1[[j]],Λ]-CWeyl[3,λs[[j]],1]Vol3[1]  Λ^(3/2), Bfree[3,αλ[λs[[j]],1],1]Area3[1]Λ,-BDirLiu[3,λs[[j]],1]Area3[1]Λ},{Λ, 0, Lmax1[[j]]}, PlotStyle->{Black,Blue,Red},Exclusions->None,
Axes->True,
AxesLabel->{MaTeX["\\Lambda"], None},
PlotRange->{{0,Lmax1[[j]]},Automatic},
Ticks->{xticks1[[j]],yticks1},
Epilog->Inset[Framed[MaTeX["\\lambda="<>lambdatext[[j]]<>",\\quad h=1"], RoundingRadius->2, ContentPadding->False,FrameMargins->Tiny],Scaled[{0.1,1}],Scaled[{0,1}]]
],{j,1,3}]

legendcylfreenew=LineLegend[{Black, Blue, Red},MaTeX[{"\\text{computed }\\mathcal{N}_{\\mathrm{free}}(\\Lambda)- a (2\\pi)^2 h \\Lambda^{3/2}\\qquad\\qquad","8\\pi^2 b_{\\mathrm{free}}\\Lambda","8\\pi^2 b_{\\mathrm{free}}^{\\mathrm{Liu}}\\Lambda"}],LegendFunction->None,LegendLayout->{"Row",1}]

figcylfree=GraphicsColumn[{GraphicsRow[graphicstable1Free,ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
legendcylfree},ImageSize->Full,Spacings->{0,Automatic}]

Export["figcylfree.pdf",figcylfree]

figcylfreenew=GraphicsColumn[{GraphicsRow[graphicstable1Freenew,ImageSize->Full,Spacings->{Scaled[0.05],Automatic}],
legendcylfreenew},ImageSize->Full,Spacings->{0,Automatic}]

Export["figcylfreenew.pdf",figcylfreenew]