Rz

Clear[Reom,Rtest,myR,Rzic,myseries,eqlist,myequrels,myequseries,r1,T];
r1[a_,q_]=r/.Solve[((r^4)+(a^2) ((L^2)+(r^2))+(L^2) ((q^2)-2 M r+(r^2))==0)/.{L->1,M->1},r][[4]];
z1[a_,q_]:=1/r1[a,q];
r1[.5,.6]
z1[.5,.6]
T[a_,q_]:=Simplify[r1[a,q](((1+(a^2)(L^(-2))+3((r1[a,q])^2)(L^(-2))-((a^2)+(q^2))((r1[a,q])^(-2))))/(4π((a^2)+((r1[a,q])^2))))/.L->1]

{Reom[a_,q_,λ_,Z_]=Simplify[(((4 (a^2)+4 (r^2)+λ) Rr[r])/(2 ((q^2)+(a^2) (1+(r^2))+r (-2+r+(r^3)))))+((2 (-1+r+(a^2) r+2 (r^3)) Rr'[r])/((q^2)+(a^2) (1+(r^2))+r (-2+r+(r^3))))+Rr''[r]/.Rr->Function[r,Rz[Z/r]/r]/.r->Z /z ]}//MatrixForm
Rtest[z_]:=1+ bb1(z-1)+ cc1(z-1) ^2+ dd1(z-1) ^3 + ee1(z-1) ^4;
{myseries[a_,q_,λ_,Z_]=Series[Reom[a,q,λ,Z]/.  Rz->Rtest, {z,1,3}]}//MatrixForm
(eqlist[a_,q_,λ_,Z_]=Complement[Table[Simplify[SeriesCoefficient[myseries[a,q,λ,Z],i]]==0, {i,-1,2}],{True}])//MatrixForm
(myequrels[a_,q_,λ_,Z_] = Solve[eqlist[a,q,λ,Z], {cc1,dd1,ee1}])//MatrixForm
Rzic[ a_,q_,λ_,Z_,bb1_][z_] := Evaluate[Rtest[z] /. myequrels[a,q,λ,Z][[1]]]

Clear[myfalloffs,myR]
small = 1/10000000
delta = 1/1000000
myR[a_,q_,λ_,Z_,bb1_?NumberQ] :=Block[{Rz},  Rz /. NDSolve[{Reom[a,q,λ,Z]==0,Rz[1-small]==Rzic[a,q,λ,Z,bb1][1-small],
Rz'[1-small]==Rzic[ a,q,λ,Z,bb1]'[1-small]},Rz,{z, 1- small,delta}][[1]]];
With[{a=.1,q=.5,λ=2},{Rzic[ a,q,λ,z1[a,q],406/100][1-small],myR[ a,q,λ,z1[a,q],406/100]}]

Rzzero[z_] := AA1 + BB1 z+ CC1 z^2 +DD1 z^3+EE1 z^4
{myequseries[a_,q_,λ_,Z_] = Series[Reom[a,q,λ,Z]/. Rz->Rzzero, {z,0,7}]}//MatrixForm
(eqlist[a_,q_,λ_,Z_]=Complement[Table[Simplify[SeriesCoefficient[myequseries[a,q,λ,Z],i]]==0, {i,5,7}],{True}])//MatrixForm
myequrels[a_,q_,λ_,Z_] = First[ Solve[eqlist[a,q,λ,Z], {CC1,DD1,EE1}]]

myfalloffs[a_,q_,λ_,bb1_]:={AA1,BB1}/.Solve[{myR[a,q,λ,z1[a,q],bb1][delta] == Rzzero[delta] /.myequrels[a,q,λ,z1[a,q]],
myR[a,q,λ,z1[a,q],bb1]'[delta] == Rzzero'[delta]/.myequrels[a,q,λ,z1[a,q]]},{AA1,BB1}][[1]]
With[{a=.1,q=.6,λ=2},myfalloffs[a,q,λ, 406/100]]

FindRoot[First[myfalloffs[1/100,2/100,1,yy]]==0,{yy,406/100}]

mapitone[ a_,q_,λ_, yseed_] :=Block[{yy}, Join[{q},{yy/.FindRoot[N[myfalloffs[a,q,λ,yy][[2]]]==0, {yy,yseed}]}]]
mapittwo[ a_,q_,λ_, yseed_] := Block[{yy},Join[{q},{yy/.FindRoot[N[myfalloffs[a,q,λ,yy][[1]]]==0, {yy,yseed}]}]]

With[{a=1/10,q=6/10,λ=1/10},{mapitone[a,q,λ,1],mapittwo[a,q,λ,4]}]
With[{a=1/10,q=83502078/100000000,λ=1/10},myTc1 =T[a,q];N[{T[a,q]/myTc1, Sqrt[2]myfalloffs[a,q,λ,mapitone[a,q,λ,1][[2]]]/myTc1}]]

iilist={0.00021`,0.00031`,0.00041`,0.00051`,0.00061`,0.00071`,0.00081`,0.00091`,0.00101`,0.00201`,0.00301`,0.00401`,0.00501`,0.00601`,0.00701`,0.00801`,0.00901`,0.01001`,0.02001`,0.03001`,0.04001`,0.05001`,0.06001`,0.07001`,0.08001`,0.09001`,0.10001`,0.15001`,0.20001`,0.25001`,0.30001`,0.35001`,0.40001`,0.45001`,0.50001`,0.5067004156`,0.5134008312`,0.5201012468`,0.5268016624`,0.533502078`,0.5402024936`,0.5469029092`,0.5536033248`,0.5603037404`,0.567004156`,0.5737045716`,0.5804049872`,0.5871054028`,0.5938058184`,0.600506234`,0.6072066496`,0.6139070652`,0.6206074807999999`,0.6273078964000001`,0.634008312`,0.6407087276`,0.6474091432`,0.6541095588`,0.6608099744`,0.66751039`,0.6742108056`,0.6809112212`,0.6876116368`,0.6943120523999999`,0.701012468`,0.7077128836`,0.7144132992000001`,0.7211137148`,0.7278141304000001`,0.734514546`,0.7412149616`,0.7479153772`,0.7546157928`,0.7613162084`,0.768016624`,0.7747170396`,0.7814174551999999`,0.7881178708000001`,0.7948182864`,0.8015187020000001`,0.8082191176`,0.8149195332`,0.8216199488`,0.8283203644`,0.83502078`}

Block[{a=55/1000,λ=1/10,y0=1,temprels},dimonesweep=Table[{{temprels = mapitone[a,ii,λ,y0], y0 = temprels[[2]] }[[1]],myfalloffs[a,ii,λ,temprels[[2]]]},{ii,iilist}];InputForm[dimonesweep, NumberMarks->False]]

Block[{a=1/10,λ=1/10,y0=1,temprels},dimonesweep1=Table[{{temprels = mapitone[a,ii,λ,y0], y0 = temprels[[2]] }[[1]],myfalloffs[a,ii,λ,temprels[[2]]]},{ii,iilist}];InputForm[dimonesweep1, NumberMarks->False]]

Block[{a=1/10,λ=1/10,y0=1,temprels},dimtwosweep=Table[{{temprels = mapittwo[a,ii,λ,y0],y0 = temprels[[2]]}[[1]],myfalloffs[a,ii,λ,temprels[[2]]]},{ii,iilist}];InputForm[dimtwosweep, NumberMarks->False]]

Block[{a=15/100,λ=1/10,y0=1,temprels},dimtwosweep1=Table[{{temprels = mapittwo[a,ii,λ,y0],y0 = temprels[[2]]}[[1]],myfalloffs[a,ii,λ,temprels[[2]]]},{ii,iilist}];InputForm[dimtwosweep1, NumberMarks->False]]

With[{a=55/1000},myTc1 =T[a,dimonesweep[[1]][[1]][[1]]];
p1=ListLinePlot[Table[{T[a,dimonesweep[[i]][[1]][[1]]]/myTc1, Sqrt[2]dimonesweep[[i]][[2]][[1]]/myTc1},{i,1,Length[dimonesweep]}], FrameLabel-> {T/(Subscript[T, c]), "<(Subscript[O, 1])>"/(Subscript[T, c])}, Frame->True, RotateLabel->False]]
With[{a=1/10},myTc1 =T[a,dimonesweep1[[1]][[1]][[1]]];
p2=ListLinePlot[Table[{T[a,dimonesweep1[[i]][[1]][[1]]]/myTc1, Sqrt[2]dimonesweep1[[i]][[2]][[1]]/myTc1},{i,1,Length[dimonesweep1]}], FrameLabel-> {T/(Subscript[T, c]), "<(Subscript[O, 1])>"/(Subscript[T, c])}, Frame->True, RotateLabel->False,PlotStyle->Red]]
Show[{p1,p2}]

With[{a=.1},myTc2 =T[a,dimtwosweep[[1]][[1]][[1]]];
p1=ListLinePlot[Table[{T[a,dimtwosweep[[i]][[1]][[1]]]/myTc2, Sqrt[2dimtwosweep[[i]][[2]][[2]]]/myTc2},{i,1,Length[dimtwosweep]}], FrameLabel-> {T/(Subscript[T, c]), Sqrt["<(Subscript[O, 2])>"]/(Subscript[T, c])}, Frame->True, RotateLabel->False]]
With[{a=.15},myTc2 =T[a,dimtwosweep[[1]][[1]][[1]]];
p2=ListLinePlot[Table[{T[a,dimtwosweep[[i]][[1]][[1]]]/myTc2, Sqrt[2dimtwosweep[[i]][[2]][[2]]]/myTc2},{i,1,Length[dimtwosweep]}], FrameLabel-> {T/(Subscript[T, c]), Sqrt["<(Subscript[O, 2])>"]/(Subscript[T, c])}, Frame->True, RotateLabel->False,PlotStyle->Red]]
Show[{p1,p2}]

Stheta

Clear[Seom,mySx,mys];
smallx=1/100000;
Seom[a_,λ_]=Collect[Simplify[((λ+4 (a^2) ((Sin[θ])^2)) Sθ[θ]+(-2-(a^2)+3 (a^2) Cos[2 θ]) Cot[θ] (Sθ)'[θ]+(-2+(a^2)+(a^2) Cos[2 θ]) (Sθ)''[θ])/.θ->ArcCos[x]/.Sθ->Function[θ, Sx[Cos[θ]]]],{Sx''[x],Sx'[x]}, Simplify]

With[{λ=1/10},{Seom[0,λ],DSolve[Seom[0,λ]==0,Sx,x]}]
Solve[(n^2)+n+λ==0,n]

Clear[mySx,mys]
mySx[a_,λ_,s_?NumberQ]:=Block[{Sx},Sx/.NDSolve[{Seom[a,λ]==0,Sx[0]==N[LegendreP[(1/2) (-1+(Sqrt[1-4 λ])),0]],Sx'[0]==N[s]},Sx,{x,0,1-smallx}][[1]]];
mys[a_,λ_?NumberQ]:=Block[{ss,s0,xxx},s0=N[(D[LegendreP[(1/2) (-1+(Sqrt[1-4 λ])),xxx],xxx])/.xxx->0];ss/.FindRoot[mySx[a,λ,ss][1-smallx],{ss,s0}]];
Block[{λ=1/10,x},N[D[LegendreP[(1/2) (-1+(Sqrt[1-4 λ])),x],x]/.x->0]]

Block[{a=1/10,λ=1/10,s},s=N[D[LegendreP[(1/2) (-1+(Sqrt[1-4 λ])),x],x]/.x->0];mySx[a,λ,s]]
Block[{a=1/10,λ=1/10},mys[a,λ]]
Block[{a=1/10,λ=1/10,s},s=mys[a,λ];{s,mySx[a,λ,s][1-smallx]-1}]

With[{a=0,λ=1/10},p1=Plot[mySx[a,λ,mys[a,λ]][y],{y,0,1-smallx}]]
With[{a=2/10,λ=1/10},p2=Plot[mySx[a,λ,mys[a,λ]][y],{y,0,1-smallx}]]
With[{a=4/10,λ=1/10},p3=Plot[mySx[a,λ,mys[a,λ]][y],{y,0,1-smallx}]]
With[{a=6/10,λ=1/10},p4=Plot[mySx[a,λ,mys[a,λ]][y],{y,0,1-smallx}]]
With[{a=8/10,λ=1/10},p5=Plot[mySx[a,λ,mys[a,λ]][y],{y,0,1-smallx}]]

Conductivity

Clear[myAeqseries,myAeqrels]

{Aeom[ω_,a_,q_,Z_,Rpsi_,Spsi_]=Collect[Expand[((Z^2)/(z^4))(((Ax[r,θ] ((r^2)+(a^2) ((Cos[θ])^2)) (((2 (ω^2) ((a^6)-2 (r^4)+(a^4) (-1+3 (r^2))+(a^2) ((q^2)+r (-2-3 r+2 (r^3)))+(a^2) ((a^4)-(q^2)+2 r-(r^2)+(a^2) (-1+(r^2))) Cos[2 θ]))/(((q^2)+(a^2) (1+(r^2))+r (-2+r+(r^3))) (-2+(a^2)+(a^2) Cos[2 θ]) ((a^2)+2 (r^2)+(a^2) Cos[2 θ])))-2 ((ψ[r,θ])^2)))/((q^2)-2 r+(1+(r^2)) ((a^2)+(r^2))))+((4 ((r^2)+(a^2) ((Cos[θ])^2)) ((a^4) r+(a^2) (-1-r+2 (r^3))+2 r (-(q^2)+r+(r^4))+(a^2) (-1+r+(a^2) r+2 (r^3)) Cos[2 θ]) Derivative[1,0][Ax][r,θ])/(((q^2)-2 r+(1+(r^2)) ((a^2)+(r^2))) ((((a^2)+2 (r^2)+(a^2) Cos[2 θ]))^2)))+Derivative[2,0][Ax][r,θ])]/.Ax->Function[{r,theta},Ax[Z/r,Cos[theta]]]/.θ->ArcCos[x]/.ψ->Function[{r,theta},z Rpsi Spsi]/.r->Z/z ,{Derivative[2,0][Ax][z,x],Derivative[1,0][Ax][z,x],Ax[z,x]},FullSimplify]}//MatrixForm
Atest[z_,x_]=Simplify[(1-(x^2))(((1-z))^(I ω))(1+aa(1-z)+bb(((1-z))^2)+cc(((1-z))^3))];
{myAeqseries[ ω_,a_,q_,Z_,Rpsi_,Spsi_][z_,x_]=Simplify[Series[Simplify[(((1-z))^(-I ω+1))Aeom[ω,a,q,Z,Rpsi,Spsi]/.Ax->Atest],{z,1,2}]]}//MatrixForm
(myAeqrels[ ω_,a_,q_,Z_,Rpsi_,Spsi_][z_,x_]=Simplify[Solve[Table[Simplify[SeriesCoefficient[myAeqseries[ ω,a,q,Z,Rpsi,Spsi][z,x],i],Assumptions->0<x<1]==0,{i,0,2}],{aa,bb,cc}][[1]]])//MatrixForm

Axic[ ω_,a_,q_,λ_,bb1_][z_,x_] := N[Block[{s,Rpsi,Spsi,Z},s=mys[a,λ];Rpsi=Rzic[ a,q,λ,z1[a,q],bb1][z];Spsi=mySx[a,λ,s][x];Z=z1[a,q];Atest[z,x]/.myAeqrels[ ω,a,q,Z,Rpsi,Spsi][z,x]]];
Block[{i=1,ω=10,a=55/1000,q=6/10,λ=1/10,bb1=406/100},Axic[ω,a,q,λ,bb1][1-small,0.5]]

Clear[myAx]
myAx[ω_,a_,q_,λ_,bb1_?NumberQ] := Block[{Ax},Ax /. NDSolve[{Aeom[ω,a,q,z1[a,q],myR[a,q,λ,z1[a,q],bb1][z],mySx[a,λ,mys[a,λ]][x]]==0, Ax[(1-small),x]==Axic[ ω,a,q,λ,bb1][1-small,x],Derivative[1,0][Ax][1-small,x] == Derivative[1,0][Axic[ ω,a,q,λ,bb1]][1-small,x]}, Ax, {z,(1-small), delta},{x,0,1-smallx}][[1]]]

Block[{ω=72,a=55/1000,q=dimonesweep[[1]][[1]][[1]],λ=0,bb1=dimonesweep[[1]][[1]][[2]],mA},mA=myAx[ω,a,q,λ,bb1];{Plot[Re[mA[z,0]],{z,1-small,delta}],Plot[Im[mA[z,0]],{z,1-small,delta}]}]
Block[{ω=10,a=0,q=dimonesweep[[1]][[1]][[1]],λ=0,bb1=dimonesweep[[1]][[1]][[2]],mA},mA=myAx[ω,a,q,λ,bb1];{Plot[Re[mA[z,0]],{z,1-small,delta}],Plot[Im[mA[z,0]],{z,1-small,delta}]}]

mycond[i_,ω_,xx_,a_,λ_] := Block[{q=dimonesweep[[i]][[1]][[1]],bb1=dimonesweep[[i]][[1]][[2]],mA,map,ma},mA=myAx[ω,a,q,λ,bb1];map=Derivative[1,0][mA][delta,xx];ma=mA[delta,xx];Re[I map/ma/ω]]
Block[{i=1,ω=10,xxx=0,a=55/1000,λ=1/10,mA},mycond[i,ω,xxx,a,λ]]

figure2curve1[i_,xxx_,colr_] := Block[{a=0,λ=0},ListLinePlot[Table[{ω, mycond[i,ω 3/4/Pi  ,xxx,a,λ] }, {ω,2,40, 1}], RotateLabel->False, FrameLabel->{ω/T, Re[σ]},PlotStyle->colr]]
Length[dimonesweep]

p1=figure2curve1[25,0,Green]
p2=figure2curve1[50,0,Blue]
p3=figure2curve1[85,0,Red]
figure2a = Show[{p1,p2,p3}]

p1=p3
p2=figure2curve1[85,0.25,Blue]
p3=figure2curve1[85,0.5,Green]
figure2a = Show[{p1,p2,p3}]