Get["~/doc/research/first_order_singularities/mma/schofield.wl"]

SetAttributes[untilSuccess,HoldAll];
untilSuccess[e_]:=Quiet[Check[e,untilSuccess[e]]]

randomEval[eq_]:=eq/.Around[x_,y_]:>RandomVariate[NormalDistribution[x,y]]

Matching all

n1=2;
F1=RF[n1];
hf[0]=h[0];
eqf[0]=eq[n1,hf[0]][0,0,2];

sol[0]=untilSuccess[newSol[eqf[0],{θc->1+RandomReal[],θ0->0.2,B->RandomReal[],gC[0]->1},{AL,AH,A[0],A[1],A[2]},1]]

mt=7;ListLogPlot[{
Abs[Gls],
Chop[Abs[dFdξLowList[n1,hf[0]][5]/(Range[0,5]!)]/.sol[0]],
(m|->(1/π)((OverBar[s])/(2π))(((((2(OverBar[s]))/π)))^((m-1)))Gamma[m-1])/@Range[0,12]
}]

n2=2;
F2=RF[n2];
hf[1]=h[1];
eqf[1]=eq[n2,hf[1]][1,0,2];

sol[1]=untilSuccess[newSol[eqf[1],sol[0],{gC[1]},(10^(-1))]]

n3=2;
F3=RF[n3];
hf[2]=h[2];
eqf[2]=eq[n3,hf[2]][1,1,2];

sol[2]=untilSuccess[newSol[eqf[2],sol[1],{gC[2]},(10^(-3))]]

n4=2;
F4=RF[n4];
hf[3]=h[5];
eqf[3]=eq[n4,hf[3]][3,3,4];

sol[3]=untilSuccess[newSol[eqf[3],sol[2],{gC[4],gC[5]},(10^(-3))]]

n5=2;
F5=RF[n5];
hf[4]=h[7];
eqf[4]=eq[n5,hf[4]][4,4,4];

sol[4]=untilSuccess[newSol[eqf[4],sol[3],{gC[6],gC[7]},(10^(-5))]]

n6=2;
F6=RF[n5];
hf[5]=h[9];
eqf[5]=eq[n6,hf[5]][5,4,6];

sol[5]=untilSuccess[newSol[eqf[5],sol[4],{gC[8],gC[9]},(10^(-6))]]

n7=2;
F7=RF[n7];
hf[6]=h[11];
eqf[6]=eq[n7,hf[6]][6,5,6];

sol[6]=untilSuccess[newSol[eqf[6],sol[5],{gC[10],gC[11]},(10^(-6))]]

n8=2;
F8=RF[n8];
hf[7]=h[13];
eqf[7]=eq[n8,hf[7]][6,7,6];

sol[7]=untilSuccess[newSol[eqf[7],sol[6],{gC[12],gC[13]},(10^(-7))]]

n9=2;
F9=RF[n9];
hf[8]=h[15];
eqf[8]=eq[n9,hf[8]][6,8,8];

sol[8]=untilSuccess[newSol[eqf[8],sol[7],{gC[14],gC[15]},(10^(-8))]]

n10=2;
F10=RF[n10];
hf[9]=h[17];
eqf[9]=eq[n10,hf[9]][8,8,8];

sol[9]=newSol[eqf[9],sol[8],{gC[16],gC[17]},(10^(-9)),0,WorkingPrecision->7]

coeffs=(n|->{#,Abs[SeriesCoefficient[hf[n][x θc]/.sol[n],{x,0,2#-1}]]}&/@Range[1,20])/@Range[0,8];

ListLogPlot[coeffs,Joined->True]

converge=(n|->{1/(2#-1),((SeriesCoefficient[hf[n][x θc]/.sol[n],{x,0,2#-1}])/(SeriesCoefficient[hf[n][x θc]/.sol[n],{x,0,2#-3}]))}&/@Range[2,20])/@Range[0,8];

ListPlot[converge,Joined->True,PlotRange->{-3.5,0}]

Plot[Evaluate[(hf[#][x θc]/.sol[#]&/@Range[4,8])],{x,0,1.3},PlotLegends->True]

Plot[Evaluate[((hf[#][x θc]/.sol[#]&/@Range[4,8])-hf[8][x θc]/.sol[8])/(hf[8][x θc]/.sol[8])],{x,0,1.3},PlotLegends->True]

Plot[Evaluate[F2[x θc]/.sol[#]&/@Range[4,8]],{x,0,1.2}]

(t^(2-α)) F[ξ]===(R^(2-α))((t[θ])^(2-α)) F[θ]

ParametricPlot[Evaluate[{t[x θc]((hf[#][x θc])^(-1/Δ[2])),hf[#][θ]^(-2/Δ[2]) (F2[θ]-t[θ]^2 Log[hf[#][θ]^2]/(8 π Δ[2]))}/.sol[#]&/@Range[4,8]],{x,0.01,0.99},AspectRatio->1,ImageSize->500]

dFdηList[2,hf[8]][3][1]/.sol[8]

(Φs[[1]]-hf[#][θ]^(-2/Δ[2]) (F2[θ])/.θ->1/.sol[#])&/@Range[0,8]

((dFdηList[2,hf[#]][0][1][[1]]-Φs[[1]])/.sol[#])&/@Range[0,8]

θc/.sol[#]&/@Range[4,8]

ParametricPlot[Evaluate[{-t[θ]((hf[#][θ])^(-1/Δ[2])),(ddξ[hf[#]][((hf[#][θ])^((-2/Δ[2])))(F2[θ]-t[θ]^2 Log[hf[#][θ]^2]/(8 π Δ[2]))])}/.sol[#]&/@Range[4,8]],{θ,0.001,1.2},AspectRatio->1,ImageSize->500]

ParametricPlot[Evaluate[{-t[θ]((hf[#][θ])^(-1/Δ[2])),(((hf[#][θ])^((-2/Δ[2])))(F2[θ]-t[θ]^2 Log[hf[#][θ]^2]/(8 π Δ[2]))-Evaluate[((hf[#][θ])^((-2/Δ[2])))(F2[θ]-t[θ]^2 Log[hf[#][θ]^2]/(8 π Δ[2]))/.sol[8]])}/.sol[#]&/@Range[4,8]],{θ,0.001,1.2},AspectRatio->1,ImageSize->500,PlotRange->{Automatic,{-0.02,0.01}}]

ParametricPlot[Evaluate[{hf[#][x θc]((t[x θc])^(-Δ[2])),(F2[x θc]-Evaluate[F2[x θc]/.sol[8]])}/.sol[#]&/@Range[4,8]],{x,0,1.2}]

radConvFit=LinearModelFit[{1/(#+1),Abs[((gC[#])/(gC[#-1]))]}&/@Range[3,8]/.sol[8],x,x]

Export["~/doc/research/first_order_singularities/paper/data/h_series_ours_"<>ToString[Length[#]-1]<>".dat",#]&@Select[#,MatchQ[#[[1]],gC[_]]&][[All,2]]&/@(sol/@Range[0,9])

sol[8]

Show[
Plot[radConvFit[x],{x,0,1},PlotStyle->Black],
ListPlot[{1/(#+1),Abs[((gC[#])/(gC[#-1]))]}&/@Range[1,15]/.#&/@(sol/@Range[5,8]),Joined->True]
]

ListPlot[Evaluate[θc/.sol/@Range[0,9]]]

ListLogPlot[{
Abs[Gls]((Range[Length[Gls]]-1)!),
Chop[Abs[dFdξLowList[n4,h4][9]]/.sol4]
}]

n7=2;
F7=RF[n5];
h7=h[11];
eq7=eq[n7,h7][10,10];

sol7=untilSuccess[newSol[eq7,sol6,{gC[10],gC[11]},(10^(-5))]]

n8=2;
F8=RF[n8];
h8=h[12];
eq8=eq[n8,h8][11,10];

sol8=untilSuccess[newSol[eq8,sol7,{gC[12]},(10^(-6))]]

n9=2;
F9=RF[n9];
h9=h[13];
eq9=eq[n9,h9][11,12];

sol9=untilSuccess[newSol[eq9,sol8,{gC[13]},(10^(-7))]]

Φs[[1]]

(dFdη[n1,h1][0][1]/.sol1)

sol1

nDers=3;
ListLogLogPlot[
Transpose@Abs[{
((dFdηList[n1,h1][nDers][1]/.sol1)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n2,h2][nDers][1]/.sol2)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n3,h3][nDers][1]/.sol3)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n4,h4][nDers][1]/.sol4)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n5,h5][nDers][1]/.sol5)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n6,h6][nDers][1]/.sol6)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n7,h7][nDers][1]/.sol7)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n8,h8][nDers][1]/.sol8)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]]
}],
Joined->True
]

(dFdηList[n3,h3][7][1]/.sol3)

eqMid[RF[n2],h2]/@Range[0,7]/.rules/.sol2

Limit[D[RF[2][θ],θ]/.sol1,θ->θc/.sol1]

dRFc[2][1]/.sol1

(Φs((Range[Length[Φs]]-1)!)-(dFdηList[n3,h3][7][1]/.sol3))/(Abs[Φs]((Range[Length[Φs]]-1)!))

χ1=chiSquared[2,h1][8];

χ1/.sol2

sol1x=FindMinimum[χ1,Evaluate[(sol2/.Rule->List)/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->100000]

ListLogPlot[{
Abs@Gls,
Chop[Abs@(dFdξLow[2,h[3]][#]/#!&/@Range[0,8])/.sol1x[[2]]]
}]

χ2=chiSquared[1,h[3]][8];

sol2x=FindMinimum[χ2,Evaluate[Join[sol1x[[2]]/.Rule->List,{{A[3],0},{gC[4],0},{gC[5],0}}]/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->20000]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F1][#],eqLowLHS[h1][#]}&/@Range[10]]/.rules/.sol1],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[Abs@{eqHighRHS[F1][#],eqHighLHS[h1][#]}&/@Range[11]]/.rules/.sol1],PlotRange->All]
}]

test=Integrate[((-Gls[[2]]["Value"])/(2π))x Exp[π /(2 Gls[[2]]["Value"] x)]/(x^n),{x,0,∞},Assumptions->{n>2}]

test2=Solve[0==eqLow[F1,h1]/@Range[10]/.sol1/.Gl[0]->0,Derivative[#][Gl][0]&/@Range[10]]

ruleAL[h1]/AL/.sol1

((eqMid[F1,h1][#])/(Derivative[#][F1][1]))&/@Range[0,7]/.sol1/.rules

Show[LogPlot[(1/π)test,{n,0,13}],
ListLogPlot[{Abs@Gls,({#+1,Abs[Derivative[#][Gl][0]]/#!}&/@Range[10])/.test2[[1]]}],PlotRange->All]

F2=RF[4];
h2=h[6];
eq2=eq[4,h2][8];

sol2=newSol[eq2,sol1,{A[3],A[4],gC[4],gC[5],gC[6]},.0001]

GraphicsRow[{
Plot[h1[θ]/.sol1,{θ,0,θc/.sol1}],
Plot[F1[θ]/.sol1,{θ,0,θc/.sol1}]
}]

ComplexPlot[Evaluate[RFC[2][θ]/.sol1],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol1]]

ComplexPlot[Evaluate[RFC[4][θ]/.sol2],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol2]]

((eqMid[F2,h2][#])/(Derivative[#][F2][1]))&/@Range[0,7]/.sol2/.rules

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F2][#],eqLowLHS[h2][#]}&/@Range[11]]/.rules/.sol2],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[-{eqHighRHS[F2][#],eqHighLHS[h2][#]}&/@Range[12]]/.rules/.sol2],PlotRange->All]
}]

Show[
ListLogPlot[Evaluate[Transpose[Abs@{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,6]]/.sol2]]
]

h3=h[9];
F3=RF[6];
eq3=eq[F3,h3][11];

sol3=newSol[eq3/.A[6]->0,sol2,{A[5],gC[7],gC[8],gC[9]},(10^(-5))]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.sol3,{θ,0,θc/.sol3}]
}]

Show[
ListPlot[Evaluate[Abs@Transpose[{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,5]]/.sol3],Joined->True]
]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F3][#],eqLowLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All],
ListLogPlot[Evaluate[Transpose[{eqHighRHS[F3][#],eqHighLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All]
}]

h4=h[4];
F4=RF[5];
eq4=eq[F4,h4][8];

sol4=newSol[eq4,sol3,{A[5],A[4]},0.001]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.sol4,{θ,0,θc/.sol4}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F4][#],eqLowLHS[h4][#]}&/@Range[9]]/.rules[h4]/.sol4],PlotRange->All],
ListLogPlot[Evaluate[Transpose[-{eqHighRHS[F4][#],eqHighLHS[h4][#]}&/@Range[11]]/.rules[h4]/.sol4],PlotRange->All]
}]

h5=h[5];
F5=RF[5];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[5]},0.001,0]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[t4,h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[6];
F5=RF[6];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[6],FC[6]},0.000001]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

ListLogPlot[Abs[gC/@Range[10]]/.sol5]

h6=h[6];
F6=RF[7];
eq6=eq[7,h6][10];

sol6=newSol[eq6,sol5,{FC[7],FC[6]},(10^(-4))]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

ListLogPlot[Evaluate[-Transpose[{eqLowRHS[3,h3][#],eqLowLHS[3,h3][#]}&/@Range[12]]/.rules[h3]/.sol3],PlotRange->All]

Limit[D[RF[1][θ],{θ,2}],θ->θc]//Expand

eqLowRHS[1,h1][1]

ParametricPlot[Evaluate[(R|->{R t[θ],(R^(Δ[2]))h2[θ]})/@{0.25,0.375,0.5,0.75,1}/.sol2/.θc->θc1],{θ,-θc1,θc1}]

Φη2=Table[{t[θ]((h2[θ])^(-1/Δ[2])),((h2[θ])^(-2/Δ[2]))(F2[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h2]/.sol2,{θ,0.0001,θc-0.0001/.sol2,0.0001}];

ListPlot[Φη2,Joined->True,PlotRange->{{-20,20},{-2,0}}]

Eθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]/(((R ))^1)}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,0,θc1-0.001,0.001}];

EMθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]],-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,-θc1,θc1,0.01}];

Graphics3D[{Black,Sphere[{0,0,0},0.1]}]

Show[
ListPointPlot3D[Flatten[Table[(({R (t[#1]+#3),(R^(Δ[2]))h2[#1]+(R^(1/8))#2,(R^(Δ[2]))h2[#1]-(R^(1/8))#2}/.sol2/.θc->θc1)&@@@EMθ2),{R,0,5,0.1}],1],AxesLabel->{(SubPlus[T]),(SubPlus[H]),(SubMinus[H])}]
]

ListPlot3D[Flatten[(({R t[#1]-R #3,(R^(Δ[2]))h2[#1]-(R^(Δ[2]))#2,(R^(Δ[2]))h2[#1]+(R^(Δ[2]))#2}/.sol2/.θc->θc1)&@@@EMθ2)/.R->Range[0.01,1,0.01],{1}]]

ListPlot3D[{,PlotRange->All]

Mξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[Mξ2,Joined->True]

χξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[χξ2,Joined->True]

h3=h[2];
F3=RF[2];
eq3=eq[2,h3][4];

sol3=newSol[eq3,sol2,{gC[1],gC[2],FC[2]}]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.B->ruleB[h3]/.sol3,{θ,0,θc/.sol3}]
}]

h4=h[3];
F4=RF[2];
eq4=eq[2,h4][5];

sol4=newSol[eq4,sol3,{gC[3]},1]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[4];
F5=RF[3];
eq5=eq[3,h5][6];

sol5=FindRoot[eq5,{{θc,θc/.sol4},{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],0.1RandomReal[{-1,1}]},{FC[3],0.1RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

h6=h[5];
F6=RF[3];
eq6=eq[3,h6][7];

sol6=FindRoot[eq6,{{θc,θc/.sol5},{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{FC[3],FC[3]/.sol5},{gC[5],0.1 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h6[θ]/.sol6,{θ,0,θc/.sol6}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6,{θ,0,θc/.sol6}]
}]

h7=h[6];
F7=RF[4];
eq7=eq[4,h7][8];

sol7=FindRoot[eq7,{{θc,θc/.sol6},{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{FC[3],FC[3]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],0.01 RandomReal[{-1,1}]},{FC[4],0.01 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h7[θ]/.sol7,{θ,0,θc/.sol7}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7,{θ,0,θc/.sol7}]
}]

h8=h[7];
F8=RF[4];
eq8=eq[4,h8][9];

sol8=FindRoot[eq8,{{θc,θc/.sol7},{A,A/.sol7},{gC[0],gC[0]/.sol7},{FC[0],FC[0]/.sol7},{FC[1],FC[1]/.sol7},{FC[2],FC[2]/.sol7},{gC[1],gC[1]/.sol7},{gC[2],gC[2]/.sol7},{gC[3],gC[3]/.sol7},{gC[4],gC[4]/.sol7},{FC[3],FC[3]/.sol7},{gC[5],gC[5]/.sol7},{gC[6],gC[6]/.sol7},{FC[4],FC[4]/.sol7},{gC[7],0.01RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h8[θ]/.sol8,{θ,0,θc/.sol8}],
Plot[F8[θ]/.B->ruleB[h8]/.sol8,{θ,0,θc/.sol8}]
}]

h9=h[8];
F9=RF[5];
eq9=eq[5,h9][10];

sol9=FindRoot[eq9,{{θc,θc/.sol8},{A,A/.sol8},{gC[0],gC[0]/.sol8},{FC[0],FC[0]/.sol8},{FC[1],FC[1]/.sol8},{FC[2],FC[2]/.sol8},{gC[1],gC[1]/.sol8},{gC[2],gC[2]/.sol8},{gC[3],gC[3]/.sol8},{gC[4],gC[4]/.sol8},{FC[3],FC[3]/.sol8},{gC[5],gC[5]/.sol8},{gC[6],gC[6]/.sol8},{FC[4],FC[4]/.sol8},{gC[7],gC[7]/.sol8},{gC[8],0.001 RandomReal[{-1,1}]},{FC[5],0.001 RandomReal[{-1,1}]}},MaxIterations->50000,Method->{"Newton", "StepControl" -> "TrustRegion"}]

GraphicsRow[{
Plot[h9[θ]/.sol9,{θ,0,θc/.sol9}],
Plot[F9[θ]/.B->ruleB[h8]/.sol9,{θ,0,θc/.sol9}]
}]

ListLogPlot[Abs[gC/@Range[9]]/.sol9]

Φη3=Table[{t[θ]((h3[θ])^(-1/Δ[2])),((h3[θ])^(-2/Δ[2]))(F3[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h3]/.sol3/.θc->θc1,{θ,0.0001,θc1-0.0001,0.0001}];

Mξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χ3ξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4+Δ[2]))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[{Φη2,Φη3},Joined->True,PlotRange->{{-20,20},{-2,0}}]

ListPlot[{Mξ2,Mξ3},Joined->True]

ListLogLogPlot[{χξ2,χξ3},Joined->True]

h4=h[7];
F4=RF[3];
eq4=eq[3,h4][4];

sol4=FindRoot[eq4/.θc->θc1,{{A,A/.sol3},{gC[0],gC[0]/.sol3},{FC[0],FC[0]/.sol3},{FC[1],FC[1]/.sol3},{FC[2],FC[2]/.sol3},{gC[1],gC[1]/.sol3},{gC[2],gC[2]/.sol3},{gC[3],gC[3]/.sol3},{gC[4],gC[4]/.sol3},{gC[5],gC[5]/.sol3},{gC[6],0.1RandomReal[]},{gC[7],0.1RandomReal[]},{FC[3],0.1RandomReal[]}},MaxIterations->500]

GraphicsRow[{
Plot[h4[θ]/.sol4/.θc->θc1,{θ,0,θc1}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol4]]

h5=h[8];
F5=RF[4];
eq5=eq[4,h5][5];

sol5=FindRoot[eq5/.θc->θc1,{{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],gC[4]/.sol4},{gC[5],gC[5]/.sol4},{gC[6],gC[6]/.sol4},{gC[7],gC[7]/.sol4},{FC[3],FC[3]/.sol4},{gC[8],0.1RandomReal[]},{FC[4],0.1RandomReal[]}},MaxIterations->5000]

GraphicsRow[{
Plot[h5[θ]/.sol5/.θc->θc1,{θ,0,θc1}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol5]]

h6=h[10];
F6=RF[5];
eq6=eq[5,h6][6];

sol6=FindRoot[eq6/.θc->θc1,{{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{gC[5],gC[5]/.sol5},{gC[6],gC[6]/.sol5},{gC[7],gC[7]/.sol5},{FC[3],FC[3]/.sol5},{gC[8],gC[8]/.sol5},{FC[4],FC[4]/.sol5},{gC[9],(10^(-3))RandomReal[{-1,1}]},{gC[10],(10^(-3))RandomReal[{-1,1}]},{FC[5],(10^(-3))RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h6[θ]/.sol6/.θc->θc1,{θ,0,θc1}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[10]/.sol6]]

h7=h[12];
F7=RF[5];
eq7=eq[5,h7][7];

sol7=FindRoot[eq7/.θc->θc1,{{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],gC[6]/.sol6},{gC[7],gC[7]/.sol6},{FC[3],FC[3]/.sol6},{gC[8],gC[8]/.sol6},{FC[4],FC[4]/.sol6},{gC[9],gC[9]/.sol6},{gC[10],gC[10]/.sol6},{FC[5],FC[5]/.sol6},{gC[11],0.1RandomReal[{-1,1}]},{gC[12],0.1RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h7[θ]/.sol7/.θc->θc1,{θ,0,θc1}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[14]/.sol7]]

Matching all (Subscript[θ, c]) fixed

θcf=1.1;

n1=2;
F1=RF[n1];
hf[0]=h[0];
eqf[0]=eq[n1,hf[0]][2,-1,2];

sol[0]=newSol[eqf[0]/.θc->θcf,{θ0->0.2,B->RandomReal[],gC[0]->1},{AL,AH,A[0],A[1],A[2]},1]

n2=2;
F2=RF[n2];
hf[1]=h[3];
eqf[1]=eq[n2,hf[1]][2,2,2];

sol[1]=untilSuccess[newSol[eqf[1]/.θc->θcf,sol[0],{gC[1],gC[2],gC[3]},(10^(-1))]]

n3=2;
F3=RF[n3];
hf[2]=h[5];
eqf[2]=eq[n3,hf[2]][3,3,2];

sol[2]=untilSuccess[newSol[eqf[2]/.θc->θcf,sol[1],{gC[4],gC[5]},(10^(-1))]]

n4=2;
F4=RF[n4];
hf[3]=h[7];
eqf[3]=eq[n4,hf[3]][4,3,4];

sol[3]=untilSuccess[newSol[eqf[3]/.θc->θcf,sol[2],{gC[6],gC[7]},(10^(-3))]]

n5=2;
F5=RF[n5];
hf[4]=h[9];
eqf[4]=eq[n5,hf[4]][5,4,4];

sol[4]=untilSuccess[newSol[eqf[4]/.θc->θcf,sol[3],{gC[8],gC[9]},(10^(-3))]]

n6=2;
F6=RF[n5];
hf[5]=h[9];
eqf[5]=eq[n6,hf[5]][5,4,4];

sol[5]=untilSuccess[newSol[eqf[5]/.θc->θcf,sol[4],{gC[8],gC[9]},(10^(-6))]]

n7=2;
F7=RF[n7];
hf[6]=h[11];
eqf[6]=eq[n7,hf[6]][5,5,6];

sol[6]=untilSuccess[newSol[eqf[6]/.θc->θcf,sol[5],{gC[10],gC[11]},(10^(-6))]]

n8=2;
F8=RF[n8];
hf[7]=h[13];
eqf[7]=eq[n8,hf[7]][6,6,6];

sol[7]=untilSuccess[newSol[eqf[7]/.θc->θcf,sol[6],{gC[12],gC[13]},(10^(-7))]]

n9=2;
F9=RF[n9];
hf[8]=h[15];
eqf[8]=eq[n9,hf[8]][7,7,6];

sol[8]=untilSuccess[newSol[eqf[8]/.θc->θcf,sol[7],{gC[14],gC[15]},(10^(-8))]]

n10=2;
F10=RF[n10];
hf[9]=h[17];
eqf[9]=eq[n10,hf[9]][8,8,8];

sol[9]=newSol[eqf[9],sol[8],{gC[16],gC[17]},(10^(-9)),0,WorkingPrecision->7]

test4={#,Abs[SeriesCoefficient[hf[4][θ]/.sol[4]/.θc->θcf,{θ,0,2#-1}]]}&/@Range[1,15]

test5={#,Abs[SeriesCoefficient[hf[5][θ]/.sol[5]/.θc->θcf,{θ,0,2#-1}]]}&/@Range[1,15]

test6={#,Abs[SeriesCoefficient[hf[6][θ]/.sol[6]/.θc->θcf,{θ,0,2#-1}]]}&/@Range[1,15]

test7={#,Abs[SeriesCoefficient[hf[7][θ]/.sol[7]/.θc->θcf,{θ,0,2#-1}]]}&/@Range[1,15]

test8={#,Abs[SeriesCoefficient[hf[8][θ]/.sol[8]/.θc->θcf,{θ,0,2#-1}]]}&/@Range[1,19]

ListLogPlot[{test4,test5,test6,test7,test8},Joined->True]

p[n_,k_]:=(2^n)Binomial[n,k]Binomial[(n+k-1)/2,n]

Sum[p[7,2k+1](x^(2k+1)),{k,0,4}]

p[7,1]

hk=Sum[p[2n+1,2k+1]a[n]/((θc)^k)/.rs,{n,k,∞},Assumptions->{k>1,k∈Integers,r>0,m>0}]

test=Simplify[Evaluate[hk/(hk/.k->k-1)]]/.k->1/x/.{r->0.54,m->1,θc->1.3}

Limit[Simplify[Evaluate[hn/(hn/.k->k-1)]],k->∞,Assumptions->{r>0,b>0}]

Plot[test,{x,0.005,0.5}]

(Sqrt[1./8])

1./(Sqrt[10])

LegendreP[7,x]

ftest=Sum[Evaluate[(((-1))^n)a[n]LegendreP[2n+1,x/xc]/.rs[[1]]],{n,0,200}];

Plot[ftest/.{r->(Sqrt[0.5]),m->0.675,C[1]->0.01},{x,0,1.3}]

ListPlot[a[n]/.rs[[1]]

n11=2;
F11=RF[n11];
h11=h[17];
eq11=eq[n11,h11][8,7,10];

sol11=newSol[eq11,sol10,{gC[17]},(10^(-9)),0]

Plot[Evaluate[hf[#][θ θc]/.sol[#]&/@Range[4,8]],{θ,0,1.1}]

Show[
ListPlot[Abs[((gC[#-1])/(gC[#]))]&/@Range[1,15]/.#&/@{sol1,sol2,sol3,sol4,sol5,sol6,sol7,sol8,sol9},Joined->True],
Plot[Evaluate[{θc/.sol4,θc/.sol5,θc/.sol6,θc/.sol7,θc/.sol8,θc/.sol9}],{x,0,20}]
]

radConvFit=LinearModelFit[{1/(#+1),Abs[((gC[#])/(gC[#-1]))]}&/@Range[3,8]/.sol[8],x,x]

Export["~/doc/research/first_order_singularities/paper/data/h_series_ours_"<>ToString[Length[#]-1]<>".dat",#]&@Select[#,MatchQ[#[[1]],gC[_]]&][[All,2]]&/@(sol/@Range[0,9])

sol[8]

Show[
Plot[radConvFit[x],{x,0,1},PlotStyle->Black],
ListPlot[{1/(#+1),Abs[((gC[#])/(gC[#-1]))]}&/@Range[1,15]/.#&/@(sol/@Range[5,8]),Joined->True]
]

ListPlot[Evaluate[θc/.sol/@Range[0,9]]]

ListLogPlot[{
Abs[Gls]((Range[Length[Gls]]-1)!),
Chop[Abs[dFdξLowList[n4,h4][9]]/.sol4]
}]

n7=2;
F7=RF[n5];
h7=h[11];
eq7=eq[n7,h7][10,10];

sol7=untilSuccess[newSol[eq7,sol6,{gC[10],gC[11]},(10^(-5))]]

n8=2;
F8=RF[n8];
h8=h[12];
eq8=eq[n8,h8][11,10];

sol8=untilSuccess[newSol[eq8,sol7,{gC[12]},(10^(-6))]]

n9=2;
F9=RF[n9];
h9=h[13];
eq9=eq[n9,h9][11,12];

sol9=untilSuccess[newSol[eq9,sol8,{gC[13]},(10^(-7))]]

Φs[[1]]

(dFdη[n1,h1][0][1]/.sol1)

sol1

nDers=3;
ListLogLogPlot[
Transpose@Abs[{
((dFdηList[n1,h1][nDers][1]/.sol1)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n2,h2][nDers][1]/.sol2)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n3,h3][nDers][1]/.sol3)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n4,h4][nDers][1]/.sol4)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n5,h5][nDers][1]/.sol5)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n6,h6][nDers][1]/.sol6)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n7,h7][nDers][1]/.sol7)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n8,h8][nDers][1]/.sol8)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]]
}],
Joined->True
]

(dFdηList[n3,h3][7][1]/.sol3)

eqMid[RF[n2],h2]/@Range[0,7]/.rules/.sol2

Limit[D[RF[2][θ],θ]/.sol1,θ->θc/.sol1]

dRFc[2][1]/.sol1

(Φs((Range[Length[Φs]]-1)!)-(dFdηList[n3,h3][7][1]/.sol3))/(Abs[Φs]((Range[Length[Φs]]-1)!))

χ1=chiSquared[2,h1][8];

χ1/.sol2

sol1x=FindMinimum[χ1,Evaluate[(sol2/.Rule->List)/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->100000]

ListLogPlot[{
Abs@Gls,
Chop[Abs@(dFdξLow[2,h[3]][#]/#!&/@Range[0,8])/.sol1x[[2]]]
}]

χ2=chiSquared[1,h[3]][8];

sol2x=FindMinimum[χ2,Evaluate[Join[sol1x[[2]]/.Rule->List,{{A[3],0},{gC[4],0},{gC[5],0}}]/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->20000]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F1][#],eqLowLHS[h1][#]}&/@Range[10]]/.rules/.sol1],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[Abs@{eqHighRHS[F1][#],eqHighLHS[h1][#]}&/@Range[11]]/.rules/.sol1],PlotRange->All]
}]

test=Integrate[((-Gls[[2]]["Value"])/(2π))x Exp[π /(2 Gls[[2]]["Value"] x)]/(x^n),{x,0,∞},Assumptions->{n>2}]

test2=Solve[0==eqLow[F1,h1]/@Range[10]/.sol1/.Gl[0]->0,Derivative[#][Gl][0]&/@Range[10]]

ruleAL[h1]/AL/.sol1

((eqMid[F1,h1][#])/(Derivative[#][F1][1]))&/@Range[0,7]/.sol1/.rules

Show[LogPlot[(1/π)test,{n,0,13}],
ListLogPlot[{Abs@Gls,({#+1,Abs[Derivative[#][Gl][0]]/#!}&/@Range[10])/.test2[[1]]}],PlotRange->All]

F2=RF[4];
h2=h[6];
eq2=eq[4,h2][8];

sol2=newSol[eq2,sol1,{A[3],A[4],gC[4],gC[5],gC[6]},.0001]

GraphicsRow[{
Plot[h1[θ]/.sol1,{θ,0,θc/.sol1}],
Plot[F1[θ]/.sol1,{θ,0,θc/.sol1}]
}]

ComplexPlot[Evaluate[RFC[2][θ]/.sol1],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol1]]

ComplexPlot[Evaluate[RFC[4][θ]/.sol2],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol2]]

((eqMid[F2,h2][#])/(Derivative[#][F2][1]))&/@Range[0,7]/.sol2/.rules

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F2][#],eqLowLHS[h2][#]}&/@Range[11]]/.rules/.sol2],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[-{eqHighRHS[F2][#],eqHighLHS[h2][#]}&/@Range[12]]/.rules/.sol2],PlotRange->All]
}]

Show[
ListLogPlot[Evaluate[Transpose[Abs@{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,6]]/.sol2]]
]

h3=h[9];
F3=RF[6];
eq3=eq[F3,h3][11];

sol3=newSol[eq3/.A[6]->0,sol2,{A[5],gC[7],gC[8],gC[9]},(10^(-5))]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.sol3,{θ,0,θc/.sol3}]
}]

Show[
ListPlot[Evaluate[Abs@Transpose[{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,5]]/.sol3],Joined->True]
]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F3][#],eqLowLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All],
ListLogPlot[Evaluate[Transpose[{eqHighRHS[F3][#],eqHighLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All]
}]

h4=h[4];
F4=RF[5];
eq4=eq[F4,h4][8];

sol4=newSol[eq4,sol3,{A[5],A[4]},0.001]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.sol4,{θ,0,θc/.sol4}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F4][#],eqLowLHS[h4][#]}&/@Range[9]]/.rules[h4]/.sol4],PlotRange->All],
ListLogPlot[Evaluate[Transpose[-{eqHighRHS[F4][#],eqHighLHS[h4][#]}&/@Range[11]]/.rules[h4]/.sol4],PlotRange->All]
}]

h5=h[5];
F5=RF[5];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[5]},0.001,0]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[t4,h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[6];
F5=RF[6];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[6],FC[6]},0.000001]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

ListLogPlot[Abs[gC/@Range[10]]/.sol5]

h6=h[6];
F6=RF[7];
eq6=eq[7,h6][10];

sol6=newSol[eq6,sol5,{FC[7],FC[6]},(10^(-4))]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

ListLogPlot[Evaluate[-Transpose[{eqLowRHS[3,h3][#],eqLowLHS[3,h3][#]}&/@Range[12]]/.rules[h3]/.sol3],PlotRange->All]

Limit[D[RF[1][θ],{θ,2}],θ->θc]//Expand

eqLowRHS[1,h1][1]

ParametricPlot[Evaluate[(R|->{R t[θ],(R^(Δ[2]))h2[θ]})/@{0.25,0.375,0.5,0.75,1}/.sol2/.θc->θc1],{θ,-θc1,θc1}]

Φη2=Table[{t[θ]((h2[θ])^(-1/Δ[2])),((h2[θ])^(-2/Δ[2]))(F2[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h2]/.sol2,{θ,0.0001,θc-0.0001/.sol2,0.0001}];

ListPlot[Φη2,Joined->True,PlotRange->{{-20,20},{-2,0}}]

Eθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]/(((R ))^1)}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,0,θc1-0.001,0.001}];

EMθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]],-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,-θc1,θc1,0.01}];

Graphics3D[{Black,Sphere[{0,0,0},0.1]}]

Show[
ListPointPlot3D[Flatten[Table[(({R (t[#1]+#3),(R^(Δ[2]))h2[#1]+(R^(1/8))#2,(R^(Δ[2]))h2[#1]-(R^(1/8))#2}/.sol2/.θc->θc1)&@@@EMθ2),{R,0,5,0.1}],1],AxesLabel->{(SubPlus[T]),(SubPlus[H]),(SubMinus[H])}]
]

ListPlot3D[Flatten[(({R t[#1]-R #3,(R^(Δ[2]))h2[#1]-(R^(Δ[2]))#2,(R^(Δ[2]))h2[#1]+(R^(Δ[2]))#2}/.sol2/.θc->θc1)&@@@EMθ2)/.R->Range[0.01,1,0.01],{1}]]

ListPlot3D[{,PlotRange->All]

Mξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[Mξ2,Joined->True]

χξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[χξ2,Joined->True]

h3=h[2];
F3=RF[2];
eq3=eq[2,h3][4];

sol3=newSol[eq3,sol2,{gC[1],gC[2],FC[2]}]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.B->ruleB[h3]/.sol3,{θ,0,θc/.sol3}]
}]

h4=h[3];
F4=RF[2];
eq4=eq[2,h4][5];

sol4=newSol[eq4,sol3,{gC[3]},1]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[4];
F5=RF[3];
eq5=eq[3,h5][6];

sol5=FindRoot[eq5,{{θc,θc/.sol4},{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],0.1RandomReal[{-1,1}]},{FC[3],0.1RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

h6=h[5];
F6=RF[3];
eq6=eq[3,h6][7];

sol6=FindRoot[eq6,{{θc,θc/.sol5},{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{FC[3],FC[3]/.sol5},{gC[5],0.1 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h6[θ]/.sol6,{θ,0,θc/.sol6}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6,{θ,0,θc/.sol6}]
}]

h7=h[6];
F7=RF[4];
eq7=eq[4,h7][8];

sol7=FindRoot[eq7,{{θc,θc/.sol6},{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{FC[3],FC[3]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],0.01 RandomReal[{-1,1}]},{FC[4],0.01 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h7[θ]/.sol7,{θ,0,θc/.sol7}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7,{θ,0,θc/.sol7}]
}]

h8=h[7];
F8=RF[4];
eq8=eq[4,h8][9];

sol8=FindRoot[eq8,{{θc,θc/.sol7},{A,A/.sol7},{gC[0],gC[0]/.sol7},{FC[0],FC[0]/.sol7},{FC[1],FC[1]/.sol7},{FC[2],FC[2]/.sol7},{gC[1],gC[1]/.sol7},{gC[2],gC[2]/.sol7},{gC[3],gC[3]/.sol7},{gC[4],gC[4]/.sol7},{FC[3],FC[3]/.sol7},{gC[5],gC[5]/.sol7},{gC[6],gC[6]/.sol7},{FC[4],FC[4]/.sol7},{gC[7],0.01RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h8[θ]/.sol8,{θ,0,θc/.sol8}],
Plot[F8[θ]/.B->ruleB[h8]/.sol8,{θ,0,θc/.sol8}]
}]

h9=h[8];
F9=RF[5];
eq9=eq[5,h9][10];

sol9=FindRoot[eq9,{{θc,θc/.sol8},{A,A/.sol8},{gC[0],gC[0]/.sol8},{FC[0],FC[0]/.sol8},{FC[1],FC[1]/.sol8},{FC[2],FC[2]/.sol8},{gC[1],gC[1]/.sol8},{gC[2],gC[2]/.sol8},{gC[3],gC[3]/.sol8},{gC[4],gC[4]/.sol8},{FC[3],FC[3]/.sol8},{gC[5],gC[5]/.sol8},{gC[6],gC[6]/.sol8},{FC[4],FC[4]/.sol8},{gC[7],gC[7]/.sol8},{gC[8],0.001 RandomReal[{-1,1}]},{FC[5],0.001 RandomReal[{-1,1}]}},MaxIterations->50000,Method->{"Newton", "StepControl" -> "TrustRegion"}]

GraphicsRow[{
Plot[h9[θ]/.sol9,{θ,0,θc/.sol9}],
Plot[F9[θ]/.B->ruleB[h8]/.sol9,{θ,0,θc/.sol9}]
}]

ListLogPlot[Abs[gC/@Range[9]]/.sol9]

Φη3=Table[{t[θ]((h3[θ])^(-1/Δ[2])),((h3[θ])^(-2/Δ[2]))(F3[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h3]/.sol3/.θc->θc1,{θ,0.0001,θc1-0.0001,0.0001}];

Mξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χ3ξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4+Δ[2]))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[{Φη2,Φη3},Joined->True,PlotRange->{{-20,20},{-2,0}}]

ListPlot[{Mξ2,Mξ3},Joined->True]

ListLogLogPlot[{χξ2,χξ3},Joined->True]

h4=h[7];
F4=RF[3];
eq4=eq[3,h4][4];

sol4=FindRoot[eq4/.θc->θc1,{{A,A/.sol3},{gC[0],gC[0]/.sol3},{FC[0],FC[0]/.sol3},{FC[1],FC[1]/.sol3},{FC[2],FC[2]/.sol3},{gC[1],gC[1]/.sol3},{gC[2],gC[2]/.sol3},{gC[3],gC[3]/.sol3},{gC[4],gC[4]/.sol3},{gC[5],gC[5]/.sol3},{gC[6],0.1RandomReal[]},{gC[7],0.1RandomReal[]},{FC[3],0.1RandomReal[]}},MaxIterations->500]

GraphicsRow[{
Plot[h4[θ]/.sol4/.θc->θc1,{θ,0,θc1}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol4]]

h5=h[8];
F5=RF[4];
eq5=eq[4,h5][5];

sol5=FindRoot[eq5/.θc->θc1,{{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],gC[4]/.sol4},{gC[5],gC[5]/.sol4},{gC[6],gC[6]/.sol4},{gC[7],gC[7]/.sol4},{FC[3],FC[3]/.sol4},{gC[8],0.1RandomReal[]},{FC[4],0.1RandomReal[]}},MaxIterations->5000]

GraphicsRow[{
Plot[h5[θ]/.sol5/.θc->θc1,{θ,0,θc1}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol5]]

h6=h[10];
F6=RF[5];
eq6=eq[5,h6][6];

sol6=FindRoot[eq6/.θc->θc1,{{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{gC[5],gC[5]/.sol5},{gC[6],gC[6]/.sol5},{gC[7],gC[7]/.sol5},{FC[3],FC[3]/.sol5},{gC[8],gC[8]/.sol5},{FC[4],FC[4]/.sol5},{gC[9],(10^(-3))RandomReal[{-1,1}]},{gC[10],(10^(-3))RandomReal[{-1,1}]},{FC[5],(10^(-3))RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h6[θ]/.sol6/.θc->θc1,{θ,0,θc1}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[10]/.sol6]]

h7=h[12];
F7=RF[5];
eq7=eq[5,h7][7];

sol7=FindRoot[eq7/.θc->θc1,{{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],gC[6]/.sol6},{gC[7],gC[7]/.sol6},{FC[3],FC[3]/.sol6},{gC[8],gC[8]/.sol6},{FC[4],FC[4]/.sol6},{gC[9],gC[9]/.sol6},{gC[10],gC[10]/.sol6},{FC[5],FC[5]/.sol6},{gC[11],0.1RandomReal[{-1,1}]},{gC[12],0.1RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h7[θ]/.sol7/.θc->θc1,{θ,0,θc1}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[14]/.sol7]]

Match ends

nAs=2;

s0=#1->Re[#2]&@@@untilSuccess[newSol[randomEval[eqAround[nAs,h[0]][1,-1,2]],{θc->1+RandomReal[],θ0->0.2},{AH,A[1],A[2]},.1]]

Timing[dFdξLowList[nAs,h[0]][6]/(Range[0,6]!)//.rules[h[0]]/.s0]

ListLogPlot[{
Abs[Gls],
Chop[Abs[dFdξLowList[nAs,h[0]][6]]//.rules[h[0]]/.s0]
}]

ListLogPlot[{
Abs[Ghs],
Chop[Abs@dFdξHighList[nAs,h[0]][6]//.rules[h[0]]/.s0]
}]

m0=levenburgMarquardt[res[2,h[0],(10^(-6))][2]//.rules[h[0]],s0]

m01=levenburgMarquardt[res[2,h[0],(10^(-15))][2]//.rules[h[0]],m0[[2]]]

m1=levenburgMarquardt[res[2,h[1],(10^(-7))][3]//.rules[h[1]],Append[m0[[2]],gC[1]->0]]

m11=levenburgMarquardt[res[2,h[1],(10^(-15))][3]//.rules[h[1]],m1[[2]]]

m2=levenburgMarquardt[res[2,h[3],(10^(-6))][4]//.rules[h[3]],Join[m1[[2]],{gC[2]->0,gC[3]->0}]]

m21=levenburgMarquardt[res[2,h[3],(10^(-15))][4]//.rules[h[3]],Join[m2[[2]]]]

m3=levenburgMarquardt[res[2,h[4],(10^(-5))][5]//.rules[h[4]],Join[m2[[2]],{gC[4]->0}]]

m31=levenburgMarquardt[res[2,h[4],(10^(-7))][5]//.rules[h[4]],Join[m3[[2]]]]

m32=levenburgMarquardt[res[2,h[4],(10^(-9))][5]//.rules[h[4]],Join[m31[[2]]]]

m33=levenburgMarquardt[res[2,h[4],(10^(-15))][5]//.rules[h[4]],Join[m32[[2]]]]

m4=levenburgMarquardt[res[2,h[6],(10^(-5))][6]//.rules[h[6]],Join[m31[[2]],{gC[5]->0,gC[6]->0}]]

m41=levenburgMarquardt[res[2,h[6],(10^(-5))][6]//.rules[h[6]],Join[m4[[2]],{}]]

ntest=3;n1=2;
ourΦ1=dFdηList[n1,h[0]][ntest][1]/Range[0,ntest]!//.rules[h[0]]/.m01[[2]];
ourΦ2=dFdηList[n1,h[1]][ntest][1]/Range[0,ntest]!//.rules[h[1]]/.m11[[2]];
ourΦ3=dFdηList[n1,h[3]][ntest][1]/Range[0,ntest]!//.rules[h[3]]/.m21[[2]];
ourΦ4=dFdηList[n1,h[4]][ntest][1]/Range[0,ntest]!//.rules[h[4]]/.m33[[2]];
test=(Φs[[;;ntest+1]]-#)/Φs[[;;ntest+1]]&/@{ourΦ1,ourΦ2,ourΦ3,ourΦ4};

ListLogPlot[Abs@Re@Transpose[test],Joined->True]

min0=FindMinimum[Evaluate[chiSquared[2,h[0],(10^(-6))][2]//.rules[h[0]]],s0/.Rule->List,Method->"LevenbergMarquardt"]

min1=FindMinimum[Evaluate[chiSquared[2,h[0],(10^(-15))][2]//.rules[h[0]]],min0[[2]]/.Rule->List]

sol[2]=min1[[2]];hf[2]=h[0];

min3=FindMinimum[Evaluate[chiSquared[2,h[0],(10^(-6))][3]//.rules[h[0]]],Join[{#1,Re[#2]}&@@@min0[[2]],{}],MaxIterations->1000]

min4=FindMinimum[Evaluate[chiSquared[2,h[1],(10^(-6))][3]//.rules[h[1]]],Join[{#1,Re[#2]}&@@@min3[[2]],{{gC[1],0}}],MaxIterations->(10^5)]

min5=FindMinimum[Evaluate[chiSquared[2,h[1],(10^(-15))][3]//.rules[h[1]]],Join[{#1,Re[#2]}&@@@min4[[2]],{}],MaxIterations->(10^5)]

sol[3]=min5[[2]];hf[3]=h[1];

min6=FindMinimum[chiSquared[2,h[1],(10^(-6))][4]//.rules[h[1]],Join[{#1,Re[#2]}&@@@min4[[2]],{}],MaxIterations->(10^5)]

min7=FindMinimum[chiSquared[2,h[2],(10^(-5))][4]//.rules[h[2]],Join[{#1,Re[#2]}&@@@min6[[2]],{{gC[2],0}}],MaxIterations->(10^5)]

min8=FindMinimum[chiSquared[2,h[3],(10^(-5))][4]//.rules[h[3]],Join[{#1,Re[#2]}&@@@min7[[2]],{{gC[3],0}}],MaxIterations->(10^5)]

min9=FindMinimum[chiSquared[2,h[3],(10^(-15))][4]//.rules[h[3]],Join[{#1,Re[#2]}&@@@min8[[2]],{}],MaxIterations->(10^5)]

sol[4]=min9[[2]];hf[4]=h[3];

min10=FindMinimum[chiSquared[2,h[3],(10^(-5))][5]//.rules[h[3]],Join[{#1,Re[#2]}&@@@min8[[2]],{}],MaxIterations->(10^5)]

min11=FindMinimum[chiSquared[2,h[4],(10^(-5))][5]//.rules[h[4]],Join[{#1,Re[#2]}&@@@min10[[2]],{{gC[4],0}}],MaxIterations->(10^5)]

min12=FindMinimum[chiSquared[2,h[5],(10^(-5))][5]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min11[[2]],{{gC[5],0}}],MaxIterations->(10^5)]

min13=FindMinimum[chiSquared[2,h[5],(10^(-15))][5]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min12[[2]],{}],MaxIterations->(10^5)]

sol[5]=min13[[2]];hf[5]=h[5];

min14=FindMinimum[chiSquared[2,h[5],(10^(-5))][6]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min12[[2]],{}],MaxIterations->(10^5)]

min16=FindMinimum[chiSquared[2,h[6],(10^(-15))][6]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min15[[2]],{}],MaxIterations->(10^5)]

min15=FindMinimum[chiSquared[2,h[6],(10^(-5))][6]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min14[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

min16=FindMinimum[chiSquared[2,h[7],(10^(-5))][6]//.rules[h[7]],Join[{#1,Re[#2]}&@@@min42[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

sol[6]=min15[[2]];hf[6]=h[6];

h5=h[5];

ListLogPlot[{
Abs[Gls[[;;8]]-dFdξLowList[2,h5][7]/(Range[0,7]!)]//.rules[h5]/.min42[[2]]
}]

h6=h[6];

min51=FindMinimum[chiSquared[2,h[5],(10^(-5))][7]//.rules[h[5]],Join[{#1,Re[#2]}&@@@min41[[2]],{}],MaxIterations->(10^5)]

min52=FindMinimum[chiSquared[2,h[6],(10^(-5))][7]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min51[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

min53=FindMinimum[chiSquared[2,h[7],(10^(-5))][7]//.rules[h[7]],Join[{#1,Re[#2]}&@@@min52[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

min54=FindMinimum[chiSquared[2,h[6],(10^(-4))][7]//.rules[h[6]],Join[{#1,Re[#2]}&@@@min51[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

h7=h[7];

min53=FindMinimum[chiSquared[2,h7,(10^(-13))][7]//.rules[h7],Join[{#1,Re[#2]}&@@@min52[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

min52=FindMinimum[chiSquared[2,h6,(10^(-11))][7]//.rules[h6],Join[{#1,Re[#2]}&@@@min51[[2]],{}],MaxIterations->(10^5)]

h6=h[6];

min61=FindMinimum[chiSquared[2,h6,(10^(-11))][6]//.rules[h6],Join[{#1,Re[#2]}&@@@min51[[2]],{{gC[6],0}}],MaxIterations->(10^5)]

h7=h[7];

min71=FindMinimum[chiSquared[2,h7,(10^(-11))][6]//.rules[h7],Join[{#1,Re[#2]}&@@@min61[[2]],{{gC[7],0}}],MaxIterations->(10^5)]

h1'[0]

sol1s=Table[untilSuccess[newSol[randomEval[eq1],sol1,{},0,(10^(-3))]],{200}];

n2=2;
F2=RF[n2];
h2=h[2];
eq2=eqAround[n2,h2][3,-1,2];

sol2=untilSuccess[newSol[randomEval[eq2],sol1,{gC[2]},(10^(-1)),(10^(-2))]]

sol2s=Table[untilSuccess[newSol[randomEval[eq2],sol2,{},0,(10^(-3))]],{100}];

mt=8;ListLogPlot[{
Abs[Gls],
Chop[Abs[dFdξLowList[n2,h2][6]/(Range[0,6]!)]//.rules[h2]/.sol2]
}]

n3=2;
F3=RF[n3];
h3=h[4];
eq3=eqAround[n3,h3][4,-1,4];

sol3=untilSuccess[newSol[randomEval[eq3],sol2,{gC[3],gC[4]},(10^(-2))]]

sol3s=Table[untilSuccess[newSol[randomEval[eq3],sol3,{},0,(10^(-4))]],{100}];

mt=8;ListLogPlot[{
Abs[Gls],
Chop[Abs[dFdξLowList[n3,h3][6]/(Range[0,6]!)]//.rules[h3]/.sol3]
}]

ListLogPlot[{
Abs[Ghs],
Chop[Abs[dFdξHighList[n3,h3][8]/(Range[0,8]!)]//.rules[h3]/.sol3]
}]

h4=h[5];
eq4=eqAround[n3,h4][5,-1,4];

sol4=untilSuccess[newSol[randomEval[eq4],sol3,{gC[5]},(10^(-4))]]

sol4s=Table[untilSuccess[newSol[randomEval[eq4],sol4,{},0,(10^(-5))]],{100}];

h5=h[7];
eq5=eqAround[2,h5][6,-1,6];

sol5=untilSuccess[newSol[randomEval[eq5],sol4,{gC[6],gC[7]},(10^(-4))]]

sol5s=Table[untilSuccess[newSol[randomEval[eq5],sol5,{},0,(10^(-5))]],{100}];

h6=h[8];
eq6=eqAround[n3,h6][7,-1,6];

sol6=untilSuccess[newSol[randomEval[eq6],sol5,{gC[8]},(10^(-4))]]

sol6s=Table[untilSuccess[newSol[randomEval[eq6],sol6,{},0,(10^(-4))]],{100}];

h7=h[11];
eq7=eqAround[n3,h7][8,-1,8];

sol7=untilSuccess[newSol[randomEval[eq7],sol6,{gC[10],gC[11]},(10^(-4))]]

sol7s=Table[untilSuccess[newSol[randomEval[eq7],sol7,{},0,(10^(-5))]],{100}];

Plot[Evaluate[hf[#][x θc]/.rules[hf[#]]/.sol[#]&/@Range[2,5]],{x,0,1.1}]

ntest=3;n1=2;
ourΦ1=dFdηList[n1,h[0]][ntest][1]/Range[0,ntest]!//.rules[h[0]]/.m02[[2]];
ourΦ2=dFdηList[n1,h[1]][ntest][1]/Range[0,ntest]!//.rules[h[1]]/.m11;
ourΦ3=dFdηList[n1,h[3]][ntest][1]/Range[0,ntest]!//.rules[h[3]]/.m21;
ourΦ4=dFdηList[n1,h[4]][ntest][1]/Range[0,ntest]!//.rules[h[4]]/.m31;

test=(Φs[[;;ntest+1]]-#)/Φs[[;;ntest+1]]&/@{ourΦ1,ourΦ2,ourΦ3,ourΦ4};

ntest=3;
test=(Φs[[;;ntest+1]]-dFdηList[nAs,hf[#]][ntest][1]/Range[0,ntest]!)/Φs[[;;ntest+1]]//.rules[hf[#]]/.sol[#]&/@Range[2,5];

ListLogPlot[Abs@Re@Transpose[test],Joined->True]

constrainedLevenburgMarquardt[𝒞_,cons_,vars0_,γ0_:1,δ0_:1]:=Module[
{
λ,λs,𝒞cons,g,H,vars,newVars,M,δ,γ,I,nc,nv,oldCost,newCost,oldGrad,δvars,newGrad
},
nv=Length[vars0];
nc=Length[cons];
γ=γ0;δ=δ0;
λs=Thread[Rule[Array[λ,nc],ConstantArray[0,nc]]];
vars=Join[vars0,λs];
𝒞cons=𝒞+cons.(First/@λs);
Print["Compiling the gradient..."];
g=Compile[{{x,_Real,1}},Evaluate[D[𝒞cons,{First/@vars}]/.Thread[Rule[First/@vars,Part[x,#]&/@Range[nc+nv]]]]];
Print["Compiling the Hessian..."];
H=Compile[{{x,_Real,1}},Evaluate[D[𝒞cons,{First/@vars,2}]/.Thread[Rule[First/@vars,Part[x,#]&/@Range[nc+nv]]]]];
I=ArrayFlatten[{{IdentityMatrix[nv],ConstantArray[0,{nv,nc}]},{ConstantArray[0,{nc,nv}],ConstantArray[0,{nc,nc}]}}];
Print["Beginning the algorithm."];
oldCost=(𝒞/.vars);
oldGrad=Norm[g@vars[[All,2]]];
δvars=LinearSolve[Re[H@vars[[All,2]]]+γ I,Re[g@vars[[All,2]]]];
While[oldGrad>(10^(-8)),
newVars=vars;
newVars[[All,2]]-=δ δvars;
newCost=(𝒞/.newVars);
newGrad=Norm[(g@newVars[[All,2]])];
If[newGrad<oldGrad,
γ/=1.1;δ*=2;vars=newVars;oldCost=newCost;
oldGrad=newGrad;δvars=LinearSolve[Re[H@vars[[All,2]]]+γ I,Re[g@vars[[All,2]]]];,
δ/=2;];
Print[δ," ",γ," ",oldGrad," ",oldCost];
];
vars
]

n7=2;
F7=RF[n7];
h7=h[11];
eq7=eq[n7,h7][10,-1,10];

sol7=untilSuccess[newSol[eq7,sol6,{gC[10],gC[11]},(10^(-5))]]

n8=2;
F8=RF[n8];
h8=h[13];
eq8=eq[n8,h8][6,7,6];

sol8=untilSuccess[newSol[eq8,sol7,{gC[12],gC[13]},(10^(-7))]]

eqLowRHSReg[n8][#]/#!&/@Range[0,8] /.sol8

n9=2;
F9=RF[n9];
h9=h[15];
eq9=eq[n9,h9][6,8,8];

sol9=untilSuccess[newSol[eq9,sol8,{gC[14],gC[15]},(10^(-8))]]

n10=2;
F10=RF[n10];
h10=h[17];
eq10=eq[n10,h10][7,9,8];

sol10=newSol[eq10,sol9,{gC[16],gC[17]},(10^(-9))]

;;

sol10={θc->1.3935132806446657`,θ0->0.1979360333172594`,B->1.6231566110792626`,gC[0]->1.23799641800309`,AL->-0.015199741081507956`,AH->-2.470040282958811`,A[0]->1.4232276900074752`,A[1]->2.6988014458471197`,A[2]->0.24593481085642427`,gC[1]->-0.07871189880186376`,gC[2]->0.012851221284122609`,gC[3]->-0.003661937147247971`,gC[4]->0.0012371573784660372`,gC[5]->-0.00046768962383645865`,gC[6]->0.00018976673102660217`,gC[7]->-0.0000810879259692126`,gC[8]->0.00003575781335765651`,gC[9]->-0.000015998395864839337`,gC[10]->7.078544247453995`*^-6,gC[11]->-3.0030648761920763`*^-6,gC[12]->1.1728258871990909`*^-6,gC[13]->-4.014183911990179`*^-7,gC[14]->1.1268375662816859`*^-7,gC[15]->-2.3009766098169766`*^-8,gC[16]->2.5115384963668678`*^-9};

n11=2;
F11=RF[n11];
h11=h[17];
eq11=eq[n11,h11][8,7,10];

sol11=newSol[eq11,sol10,{gC[17]},(10^(-9)),0]

Plot[Evaluate[{h6[θ]/.sol6,h7[θ]/.sol7,h8[θ]/.sol8,h9[θ]/.sol9,h10[θ]/.sol10}],{θ,0,3},PlotRange->{-1,1}]

Show[
ListPlot[Abs[((gC[#-1])/(gC[#]))]&/@Range[1,15]/.#&/@{sol1,sol2,sol3,sol4,sol5,sol6,sol7,sol8,sol9},Joined->True],
Plot[Evaluate[{θc/.sol4,θc/.sol5,θc/.sol6,θc/.sol7,θc/.sol8,θc/.sol9}],{x,0,20}]
]

radConvFit=LinearModelFit[{1/#,Abs[((gC[#])/(gC[#-1]))]}&/@Range[3,8]/.sol9,x,x]

Show[
Plot[radConvFit[x],{x,0,1},PlotStyle->Black],
ListPlot[{1/#,Abs[((gC[#])/(gC[#-1]))]}&/@Range[1,15]/.#&/@{sol6,sol7,sol8,sol9},Joined->True]
]

ListPlot[Evaluate[{θc/.sol4,θc/.sol5,θc/.sol6,θc/.sol7,θc/.sol8,θc/.sol9}]]

ListLogPlot[{
Abs[Gls]((Range[Length[Gls]]-1)!),
Chop[Abs[dFdξLowList[n4,h4][9]]/.sol4]
}]

n7=2;
F7=RF[n5];
h7=h[11];
eq7=eq[n7,h7][10,10];

sol7=untilSuccess[newSol[eq7,sol6,{gC[10],gC[11]},(10^(-5))]]

n8=2;
F8=RF[n8];
h8=h[12];
eq8=eq[n8,h8][11,10];

sol8=untilSuccess[newSol[eq8,sol7,{gC[12]},(10^(-6))]]

n9=2;
F9=RF[n9];
h9=h[13];
eq9=eq[n9,h9][11,12];

sol9=untilSuccess[newSol[eq9,sol8,{gC[13]},(10^(-7))]]

Φs[[1]]

(dFdη[n1,h1][0][1]/.sol1)

sol1

nDers=3;
ListLogLogPlot[
Transpose@Abs[{
((dFdηList[n1,h1][nDers][1]/.sol1)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n2,h2][nDers][1]/.sol2)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n3,h3][nDers][1]/.sol3)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n4,h4][nDers][1]/.sol4)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n5,h5][nDers][1]/.sol5)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n6,h6][nDers][1]/.sol6)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n7,h7][nDers][1]/.sol7)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n8,h8][nDers][1]/.sol8)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]]
}],
Joined->True
]

(dFdηList[n3,h3][7][1]/.sol3)

eqMid[RF[n2],h2]/@Range[0,7]/.rules/.sol2

Limit[D[RF[2][θ],θ]/.sol1,θ->θc/.sol1]

dRFc[2][1]/.sol1

(Φs((Range[Length[Φs]]-1)!)-(dFdηList[n3,h3][7][1]/.sol3))/(Abs[Φs]((Range[Length[Φs]]-1)!))

χ1=chiSquared[2,h1][8];

χ1/.sol2

sol1x=FindMinimum[χ1,Evaluate[(sol2/.Rule->List)/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->100000]

ListLogPlot[{
Abs@Gls,
Chop[Abs@(dFdξLow[2,h[3]][#]/#!&/@Range[0,8])/.sol1x[[2]]]
}]

χ2=chiSquared[1,h[3]][8];

sol2x=FindMinimum[χ2,Evaluate[Join[sol1x[[2]]/.Rule->List,{{A[3],0},{gC[4],0},{gC[5],0}}]/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->20000]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F1][#],eqLowLHS[h1][#]}&/@Range[10]]/.rules/.sol1],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[Abs@{eqHighRHS[F1][#],eqHighLHS[h1][#]}&/@Range[11]]/.rules/.sol1],PlotRange->All]
}]

test=Integrate[((-Gls[[2]]["Value"])/(2π))x Exp[π /(2 Gls[[2]]["Value"] x)]/(x^n),{x,0,∞},Assumptions->{n>2}]

test2=Solve[0==eqLow[F1,h1]/@Range[10]/.sol1/.Gl[0]->0,Derivative[#][Gl][0]&/@Range[10]]

ruleAL[h1]/AL/.sol1

((eqMid[F1,h1][#])/(Derivative[#][F1][1]))&/@Range[0,7]/.sol1/.rules

Show[LogPlot[(1/π)test,{n,0,13}],
ListLogPlot[{Abs@Gls,({#+1,Abs[Derivative[#][Gl][0]]/#!}&/@Range[10])/.test2[[1]]}],PlotRange->All]

F2=RF[4];
h2=h[6];
eq2=eq[4,h2][8];

sol2=newSol[eq2,sol1,{A[3],A[4],gC[4],gC[5],gC[6]},.0001]

GraphicsRow[{
Plot[h1[θ]/.sol1,{θ,0,θc/.sol1}],
Plot[F1[θ]/.sol1,{θ,0,θc/.sol1}]
}]

ComplexPlot[Evaluate[RFC[2][θ]/.sol1],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol1]]

ComplexPlot[Evaluate[RFC[4][θ]/.sol2],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol2]]

((eqMid[F2,h2][#])/(Derivative[#][F2][1]))&/@Range[0,7]/.sol2/.rules

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F2][#],eqLowLHS[h2][#]}&/@Range[11]]/.rules/.sol2],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[-{eqHighRHS[F2][#],eqHighLHS[h2][#]}&/@Range[12]]/.rules/.sol2],PlotRange->All]
}]

Show[
ListLogPlot[Evaluate[Transpose[Abs@{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,6]]/.sol2]]
]

h3=h[9];
F3=RF[6];
eq3=eq[F3,h3][11];

sol3=newSol[eq3/.A[6]->0,sol2,{A[5],gC[7],gC[8],gC[9]},(10^(-5))]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.sol3,{θ,0,θc/.sol3}]
}]

Show[
ListPlot[Evaluate[Abs@Transpose[{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,5]]/.sol3],Joined->True]
]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F3][#],eqLowLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All],
ListLogPlot[Evaluate[Transpose[{eqHighRHS[F3][#],eqHighLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All]
}]

h4=h[4];
F4=RF[5];
eq4=eq[F4,h4][8];

sol4=newSol[eq4,sol3,{A[5],A[4]},0.001]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.sol4,{θ,0,θc/.sol4}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F4][#],eqLowLHS[h4][#]}&/@Range[9]]/.rules[h4]/.sol4],PlotRange->All],
ListLogPlot[Evaluate[Transpose[-{eqHighRHS[F4][#],eqHighLHS[h4][#]}&/@Range[11]]/.rules[h4]/.sol4],PlotRange->All]
}]

h5=h[5];
F5=RF[5];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[5]},0.001,0]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[t4,h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[6];
F5=RF[6];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[6],FC[6]},0.000001]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

ListLogPlot[Abs[gC/@Range[10]]/.sol5]

h6=h[6];
F6=RF[7];
eq6=eq[7,h6][10];

sol6=newSol[eq6,sol5,{FC[7],FC[6]},(10^(-4))]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

ListLogPlot[Evaluate[-Transpose[{eqLowRHS[3,h3][#],eqLowLHS[3,h3][#]}&/@Range[12]]/.rules[h3]/.sol3],PlotRange->All]

Limit[D[RF[1][θ],{θ,2}],θ->θc]//Expand

eqLowRHS[1,h1][1]

ParametricPlot[Evaluate[(R|->{R t[θ],(R^(Δ[2]))h2[θ]})/@{0.25,0.375,0.5,0.75,1}/.sol2/.θc->θc1],{θ,-θc1,θc1}]

Φη2=Table[{t[θ]((h2[θ])^(-1/Δ[2])),((h2[θ])^(-2/Δ[2]))(F2[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h2]/.sol2,{θ,0.0001,θc-0.0001/.sol2,0.0001}];

ListPlot[Φη2,Joined->True,PlotRange->{{-20,20},{-2,0}}]

Eθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]/(((R ))^1)}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,0,θc1-0.001,0.001}];

EMθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]],-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,-θc1,θc1,0.01}];

Graphics3D[{Black,Sphere[{0,0,0},0.1]}]

Show[
ListPointPlot3D[Flatten[Table[(({R (t[#1]+#3),(R^(Δ[2]))h2[#1]+(R^(1/8))#2,(R^(Δ[2]))h2[#1]-(R^(1/8))#2}/.sol2/.θc->θc1)&@@@EMθ2),{R,0,5,0.1}],1],AxesLabel->{(SubPlus[T]),(SubPlus[H]),(SubMinus[H])}]
]

Match ends (Subscript[θ, c]) fixed

θcf=1.1;

n1=2;
F1=RF[n1];
h1=h[0];
eq1=eq[n1,h1][2,-1,2];

sol1=untilSuccess[newSol[eq1/.θc->θcf,{θ0->0.2,B->RandomReal[],gC[0]->1},{AL,AH,A[0],A[1],A[2]},1]]

ListLogPlot[{
Abs[Gls]((Range[Length[Gls]]-1)!),
Chop[Abs[dFdξLowList[n1,h1][7]]/.sol1/.θc->θcf]
}]

n2=2;
F2=RF[n2];
h2=h[3];
eq2=eq[n2,h2][4,-1,4];

sol2=untilSuccess[newSol[eq2/.θc->θcf,sol1,{gC[1],gC[2],gC[3]},(10^(-1))]]

n3=2;
F3=RF[n3];
h3=h[5];
eq3=eq[n3,h3][5,-1,6];

sol3=untilSuccess[newSol[eq3/.θc->θcf,sol2,{gC[4],gC[5]},(10^(-3))]]

n4=2;
F4=RF[n4];
h4=h[7];
eq4=eq[n4,h4][6,-1,8];

sol4=untilSuccess[newSol[eq4/.θc->θcf,sol3,{gC[6],gC[7]},(10^(-3))]]

n5=2;
F5=RF[n5];
h5=h[9];
eq5=eq[n5,h5][8,-1,8];

sol5=untilSuccess[newSol[eq5,sol4,{gC[8],gC[9]},(10^(-5))]]

n6=2;
F6=RF[n5];
h6=h[9];
eq6=eq[n6,h6][9,-1,8];

sol6=untilSuccess[newSol[eq6,sol5,{gC[8],gC[9]},(10^(-6))]]

ntest=4;
ourΦ1=dFdηList[2,h1][ntest][1]/Range[0,ntest]!/.sol1;
ourΦ2=dFdηList[2,h2][ntest][1]/Range[0,ntest]!/.sol2;
ourΦ3=dFdηList[2,h3][ntest][1]/Range[0,ntest]!/.sol3;
ourΦ4=dFdηList[2,h4][ntest][1]/Range[0,ntest]!/.sol4;
ourΦ5=dFdηList[2,h5][ntest][1]/Range[0,ntest]!/.sol5;
ourΦ6=dFdηList[2,h6][ntest][1]/Range[0,ntest]!/.sol6;
ourΦ7=dFdηList[2,h7][ntest][1]/Range[0,ntest]!/.sol7;

test=(Φs[[;;ntest+1]]-#)/Φs[[;;ntest+1]]&/@{ourΦ1,ourΦ2,ourΦ3,ourΦ4,ourΦ5,ourΦ6,ourΦ7};

ListLogLogPlot[Abs@Transpose[test][[1]],Joined->True]

n7=2;
F7=RF[n7];
h7=h[11];
eq7=eq[n7,h7][10,-1,10];

sol7=untilSuccess[newSol[eq7,sol6,{gC[10],gC[11]},(10^(-5))]]

n8=2;
F8=RF[n8];
h8=h[13];
eq8=eq[n8,h8][6,7,6];

sol8=untilSuccess[newSol[eq8,sol7,{gC[12],gC[13]},(10^(-7))]]

eqLowRHSReg[n8][#]/#!&/@Range[0,8] /.sol8

n9=2;
F9=RF[n9];
h9=h[15];
eq9=eq[n9,h9][6,8,8];

sol9=untilSuccess[newSol[eq9,sol8,{gC[14],gC[15]},(10^(-8))]]

n10=2;
F10=RF[n10];
h10=h[17];
eq10=eq[n10,h10][7,9,8];

sol10=newSol[eq10,sol9,{gC[16],gC[17]},(10^(-9))]

;;

sol10={θc->1.3935132806446657`,θ0->0.1979360333172594`,B->1.6231566110792626`,gC[0]->1.23799641800309`,AL->-0.015199741081507956`,AH->-2.470040282958811`,A[0]->1.4232276900074752`,A[1]->2.6988014458471197`,A[2]->0.24593481085642427`,gC[1]->-0.07871189880186376`,gC[2]->0.012851221284122609`,gC[3]->-0.003661937147247971`,gC[4]->0.0012371573784660372`,gC[5]->-0.00046768962383645865`,gC[6]->0.00018976673102660217`,gC[7]->-0.0000810879259692126`,gC[8]->0.00003575781335765651`,gC[9]->-0.000015998395864839337`,gC[10]->7.078544247453995`*^-6,gC[11]->-3.0030648761920763`*^-6,gC[12]->1.1728258871990909`*^-6,gC[13]->-4.014183911990179`*^-7,gC[14]->1.1268375662816859`*^-7,gC[15]->-2.3009766098169766`*^-8,gC[16]->2.5115384963668678`*^-9};

n11=2;
F11=RF[n11];
h11=h[17];
eq11=eq[n11,h11][8,7,10];

sol11=newSol[eq11,sol10,{gC[17]},(10^(-9)),0]

Plot[Evaluate[{h6[θ]/.sol6,h7[θ]/.sol7,h8[θ]/.sol8,h9[θ]/.sol9,h10[θ]/.sol10}],{θ,0,3},PlotRange->{-1,1}]

Show[
ListPlot[Abs[((gC[#-1])/(gC[#]))]&/@Range[1,15]/.#&/@{sol1,sol2,sol3,sol4,sol5,sol6,sol7,sol8,sol9},Joined->True],
Plot[Evaluate[{θc/.sol4,θc/.sol5,θc/.sol6,θc/.sol7,θc/.sol8,θc/.sol9}],{x,0,20}]
]

radConvFit=LinearModelFit[{1/#,Abs[((gC[#])/(gC[#-1]))]}&/@Range[3,8]/.sol9,x,x]

Show[
Plot[radConvFit[x],{x,0,1},PlotStyle->Black],
ListPlot[{1/#,Abs[((gC[#])/(gC[#-1]))]}&/@Range[1,15]/.#&/@{sol6,sol7,sol8,sol9},Joined->True]
]

ListPlot[Evaluate[{θc/.sol4,θc/.sol5,θc/.sol6,θc/.sol7,θc/.sol8,θc/.sol9}]]

ListLogPlot[{
Abs[Gls]((Range[Length[Gls]]-1)!),
Chop[Abs[dFdξLowList[n4,h4][9]]/.sol4]
}]

n7=2;
F7=RF[n5];
h7=h[11];
eq7=eq[n7,h7][10,10];

sol7=untilSuccess[newSol[eq7,sol6,{gC[10],gC[11]},(10^(-5))]]

n8=2;
F8=RF[n8];
h8=h[12];
eq8=eq[n8,h8][11,10];

sol8=untilSuccess[newSol[eq8,sol7,{gC[12]},(10^(-6))]]

n9=2;
F9=RF[n9];
h9=h[13];
eq9=eq[n9,h9][11,12];

sol9=untilSuccess[newSol[eq9,sol8,{gC[13]},(10^(-7))]]

Φs[[1]]

(dFdη[n1,h1][0][1]/.sol1)

sol1

nDers=3;
ListLogLogPlot[
Transpose@Abs[{
((dFdηList[n1,h1][nDers][1]/.sol1)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n2,h2][nDers][1]/.sol2)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n3,h3][nDers][1]/.sol3)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n4,h4][nDers][1]/.sol4)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n5,h5][nDers][1]/.sol5)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n6,h6][nDers][1]/.sol6)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n7,h7][nDers][1]/.sol7)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]],
((dFdηList[n8,h8][nDers][1]/.sol8)/(Range[nDers+1]-1)!-Φs[[;;nDers+1]])/Φs[[;;nDers+1]]
}],
Joined->True
]

(dFdηList[n3,h3][7][1]/.sol3)

eqMid[RF[n2],h2]/@Range[0,7]/.rules/.sol2

Limit[D[RF[2][θ],θ]/.sol1,θ->θc/.sol1]

dRFc[2][1]/.sol1

(Φs((Range[Length[Φs]]-1)!)-(dFdηList[n3,h3][7][1]/.sol3))/(Abs[Φs]((Range[Length[Φs]]-1)!))

χ1=chiSquared[2,h1][8];

χ1/.sol2

sol1x=FindMinimum[χ1,Evaluate[(sol2/.Rule->List)/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->100000]

ListLogPlot[{
Abs@Gls,
Chop[Abs@(dFdξLow[2,h[3]][#]/#!&/@Range[0,8])/.sol1x[[2]]]
}]

χ2=chiSquared[1,h[3]][8];

sol2x=FindMinimum[χ2,Evaluate[Join[sol1x[[2]]/.Rule->List,{{A[3],0},{gC[4],0},{gC[5],0}}]/.{x_,y_}:>{x,y+0.000001 RandomVariate[NormalDistribution[]]}],MaxIterations->20000]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F1][#],eqLowLHS[h1][#]}&/@Range[10]]/.rules/.sol1],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[Abs@{eqHighRHS[F1][#],eqHighLHS[h1][#]}&/@Range[11]]/.rules/.sol1],PlotRange->All]
}]

test=Integrate[((-Gls[[2]]["Value"])/(2π))x Exp[π /(2 Gls[[2]]["Value"] x)]/(x^n),{x,0,∞},Assumptions->{n>2}]

test2=Solve[0==eqLow[F1,h1]/@Range[10]/.sol1/.Gl[0]->0,Derivative[#][Gl][0]&/@Range[10]]

ruleAL[h1]/AL/.sol1

((eqMid[F1,h1][#])/(Derivative[#][F1][1]))&/@Range[0,7]/.sol1/.rules

Show[LogPlot[(1/π)test,{n,0,13}],
ListLogPlot[{Abs@Gls,({#+1,Abs[Derivative[#][Gl][0]]/#!}&/@Range[10])/.test2[[1]]}],PlotRange->All]

F2=RF[4];
h2=h[6];
eq2=eq[4,h2][8];

sol2=newSol[eq2,sol1,{A[3],A[4],gC[4],gC[5],gC[6]},.0001]

GraphicsRow[{
Plot[h1[θ]/.sol1,{θ,0,θc/.sol1}],
Plot[F1[θ]/.sol1,{θ,0,θc/.sol1}]
}]

ComplexPlot[Evaluate[RFC[2][θ]/.sol1],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol1]]

ComplexPlot[Evaluate[RFC[4][θ]/.sol2],{θ,-2.5-0.5I,2.5+0.5I},AspectRatio->1/GoldenRatio,Exclusions->Evaluate[{{Im[θ]==0,((Re[θ])^2)>((θc)^2)},{Re[θ]==0,((Im[θ])^2)>((θ0)^2)}}/.sol2]]

((eqMid[F2,h2][#])/(Derivative[#][F2][1]))&/@Range[0,7]/.sol2/.rules

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F2][#],eqLowLHS[h2][#]}&/@Range[11]]/.rules/.sol2],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[-{eqHighRHS[F2][#],eqHighLHS[h2][#]}&/@Range[12]]/.rules/.sol2],PlotRange->All]
}]

Show[
ListLogPlot[Evaluate[Transpose[Abs@{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,6]]/.sol2]]
]

h3=h[9];
F3=RF[6];
eq3=eq[F3,h3][11];

sol3=newSol[eq3/.A[6]->0,sol2,{A[5],gC[7],gC[8],gC[9]},(10^(-5))]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.sol3,{θ,0,θc/.sol3}]
}]

Show[
ListPlot[Evaluate[Abs@Transpose[{gC[#]/gC[#+1],A[#]/A[#+1]}&/@Range[0,5]]/.sol3],Joined->True]
]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F3][#],eqLowLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All],
ListLogPlot[Evaluate[Transpose[{eqHighRHS[F3][#],eqHighLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All]
}]

h4=h[4];
F4=RF[5];
eq4=eq[F4,h4][8];

sol4=newSol[eq4,sol3,{A[5],A[4]},0.001]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.sol4,{θ,0,θc/.sol4}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F4][#],eqLowLHS[h4][#]}&/@Range[9]]/.rules[h4]/.sol4],PlotRange->All],
ListLogPlot[Evaluate[Transpose[-{eqHighRHS[F4][#],eqHighLHS[h4][#]}&/@Range[11]]/.rules[h4]/.sol4],PlotRange->All]
}]

h5=h[5];
F5=RF[5];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[5]},0.001,0]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[t4,h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[6];
F5=RF[6];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[6],FC[6]},0.000001]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

ListLogPlot[Abs[gC/@Range[10]]/.sol5]

h6=h[6];
F6=RF[7];
eq6=eq[7,h6][10];

sol6=newSol[eq6,sol5,{FC[7],FC[6]},(10^(-4))]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

ListLogPlot[Evaluate[-Transpose[{eqLowRHS[3,h3][#],eqLowLHS[3,h3][#]}&/@Range[12]]/.rules[h3]/.sol3],PlotRange->All]

Limit[D[RF[1][θ],{θ,2}],θ->θc]//Expand

eqLowRHS[1,h1][1]

ParametricPlot[Evaluate[(R|->{R t[θ],(R^(Δ[2]))h2[θ]})/@{0.25,0.375,0.5,0.75,1}/.sol2/.θc->θc1],{θ,-θc1,θc1}]

Φη2=Table[{t[θ]((h2[θ])^(-1/Δ[2])),((h2[θ])^(-2/Δ[2]))(F2[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h2]/.sol2,{θ,0.0001,θc-0.0001/.sol2,0.0001}];

ListPlot[Φη2,Joined->True,PlotRange->{{-20,20},{-2,0}}]

Eθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]/(((R ))^1)}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,0,θc1-0.001,0.001}];

EMθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]],-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,-θc1,θc1,0.01}];

Graphics3D[{Black,Sphere[{0,0,0},0.1]}]

Show[
ListPointPlot3D[Flatten[Table[(({R (t[#1]+#3),(R^(Δ[2]))h2[#1]+(R^(1/8))#2,(R^(Δ[2]))h2[#1]-(R^(1/8))#2}/.sol2/.θc->θc1)&@@@EMθ2),{R,0,5,0.1}],1],AxesLabel->{(SubPlus[T]),(SubPlus[H]),(SubMinus[H])}]
]

ListPlot3D[Flatten[(({R t[#1]-R #3,(R^(Δ[2]))h2[#1]-(R^(Δ[2]))#2,(R^(Δ[2]))h2[#1]+(R^(Δ[2]))#2}/.sol2/.θc->θc1)&@@@EMθ2)/.R->Range[0.01,1,0.01],{1}]]

ListPlot3D[{,PlotRange->All]

Mξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[Mξ2,Joined->True]

χξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[χξ2,Joined->True]

h3=h[2];
F3=RF[2];
eq3=eq[2,h3][4];

sol3=newSol[eq3,sol2,{gC[1],gC[2],FC[2]}]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.B->ruleB[h3]/.sol3,{θ,0,θc/.sol3}]
}]

h4=h[3];
F4=RF[2];
eq4=eq[2,h4][5];

sol4=newSol[eq4,sol3,{gC[3]},1]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[4];
F5=RF[3];
eq5=eq[3,h5][6];

sol5=FindRoot[eq5,{{θc,θc/.sol4},{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],0.1RandomReal[{-1,1}]},{FC[3],0.1RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

h6=h[5];
F6=RF[3];
eq6=eq[3,h6][7];

sol6=FindRoot[eq6,{{θc,θc/.sol5},{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{FC[3],FC[3]/.sol5},{gC[5],0.1 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h6[θ]/.sol6,{θ,0,θc/.sol6}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6,{θ,0,θc/.sol6}]
}]

h7=h[6];
F7=RF[4];
eq7=eq[4,h7][8];

sol7=FindRoot[eq7,{{θc,θc/.sol6},{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{FC[3],FC[3]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],0.01 RandomReal[{-1,1}]},{FC[4],0.01 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h7[θ]/.sol7,{θ,0,θc/.sol7}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7,{θ,0,θc/.sol7}]
}]

h8=h[7];
F8=RF[4];
eq8=eq[4,h8][9];

sol8=FindRoot[eq8,{{θc,θc/.sol7},{A,A/.sol7},{gC[0],gC[0]/.sol7},{FC[0],FC[0]/.sol7},{FC[1],FC[1]/.sol7},{FC[2],FC[2]/.sol7},{gC[1],gC[1]/.sol7},{gC[2],gC[2]/.sol7},{gC[3],gC[3]/.sol7},{gC[4],gC[4]/.sol7},{FC[3],FC[3]/.sol7},{gC[5],gC[5]/.sol7},{gC[6],gC[6]/.sol7},{FC[4],FC[4]/.sol7},{gC[7],0.01RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h8[θ]/.sol8,{θ,0,θc/.sol8}],
Plot[F8[θ]/.B->ruleB[h8]/.sol8,{θ,0,θc/.sol8}]
}]

h9=h[8];
F9=RF[5];
eq9=eq[5,h9][10];

sol9=FindRoot[eq9,{{θc,θc/.sol8},{A,A/.sol8},{gC[0],gC[0]/.sol8},{FC[0],FC[0]/.sol8},{FC[1],FC[1]/.sol8},{FC[2],FC[2]/.sol8},{gC[1],gC[1]/.sol8},{gC[2],gC[2]/.sol8},{gC[3],gC[3]/.sol8},{gC[4],gC[4]/.sol8},{FC[3],FC[3]/.sol8},{gC[5],gC[5]/.sol8},{gC[6],gC[6]/.sol8},{FC[4],FC[4]/.sol8},{gC[7],gC[7]/.sol8},{gC[8],0.001 RandomReal[{-1,1}]},{FC[5],0.001 RandomReal[{-1,1}]}},MaxIterations->50000,Method->{"Newton", "StepControl" -> "TrustRegion"}]

GraphicsRow[{
Plot[h9[θ]/.sol9,{θ,0,θc/.sol9}],
Plot[F9[θ]/.B->ruleB[h8]/.sol9,{θ,0,θc/.sol9}]
}]

ListLogPlot[Abs[gC/@Range[9]]/.sol9]

Φη3=Table[{t[θ]((h3[θ])^(-1/Δ[2])),((h3[θ])^(-2/Δ[2]))(F3[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h3]/.sol3/.θc->θc1,{θ,0.0001,θc1-0.0001,0.0001}];

Mξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χ3ξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4+Δ[2]))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[{Φη2,Φη3},Joined->True,PlotRange->{{-20,20},{-2,0}}]

ListPlot[{Mξ2,Mξ3},Joined->True]

ListLogLogPlot[{χξ2,χξ3},Joined->True]

h4=h[7];
F4=RF[3];
eq4=eq[3,h4][4];

sol4=FindRoot[eq4/.θc->θc1,{{A,A/.sol3},{gC[0],gC[0]/.sol3},{FC[0],FC[0]/.sol3},{FC[1],FC[1]/.sol3},{FC[2],FC[2]/.sol3},{gC[1],gC[1]/.sol3},{gC[2],gC[2]/.sol3},{gC[3],gC[3]/.sol3},{gC[4],gC[4]/.sol3},{gC[5],gC[5]/.sol3},{gC[6],0.1RandomReal[]},{gC[7],0.1RandomReal[]},{FC[3],0.1RandomReal[]}},MaxIterations->500]

GraphicsRow[{
Plot[h4[θ]/.sol4/.θc->θc1,{θ,0,θc1}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol4]]

h5=h[8];
F5=RF[4];
eq5=eq[4,h5][5];

sol5=FindRoot[eq5/.θc->θc1,{{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],gC[4]/.sol4},{gC[5],gC[5]/.sol4},{gC[6],gC[6]/.sol4},{gC[7],gC[7]/.sol4},{FC[3],FC[3]/.sol4},{gC[8],0.1RandomReal[]},{FC[4],0.1RandomReal[]}},MaxIterations->5000]

GraphicsRow[{
Plot[h5[θ]/.sol5/.θc->θc1,{θ,0,θc1}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol5]]

h6=h[10];
F6=RF[5];
eq6=eq[5,h6][6];

sol6=FindRoot[eq6/.θc->θc1,{{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{gC[5],gC[5]/.sol5},{gC[6],gC[6]/.sol5},{gC[7],gC[7]/.sol5},{FC[3],FC[3]/.sol5},{gC[8],gC[8]/.sol5},{FC[4],FC[4]/.sol5},{gC[9],(10^(-3))RandomReal[{-1,1}]},{gC[10],(10^(-3))RandomReal[{-1,1}]},{FC[5],(10^(-3))RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h6[θ]/.sol6/.θc->θc1,{θ,0,θc1}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[10]/.sol6]]

h7=h[12];
F7=RF[5];
eq7=eq[5,h7][7];

sol7=FindRoot[eq7/.θc->θc1,{{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],gC[6]/.sol6},{gC[7],gC[7]/.sol6},{FC[3],FC[3]/.sol6},{gC[8],gC[8]/.sol6},{FC[4],FC[4]/.sol6},{gC[9],gC[9]/.sol6},{gC[10],gC[10]/.sol6},{FC[5],FC[5]/.sol6},{gC[11],0.1RandomReal[{-1,1}]},{gC[12],0.1RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h7[θ]/.sol7/.θc->θc1,{θ,0,θc1}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[14]/.sol7]]

Caselle comparison

hCaselle[θ_]:=h0 θ(1-((θ^2)/1.16951))(1-0.222389 (θ^2)-0.043547 (θ^4)-0.014809 (θ^6)-0.007168 (θ^8))

m0Sol=Solve[-Gls[[2]]==m0 θ/((((θ^2)-1))^(β[2])) /.θ->Sqrt[1.16951],m0][[1]]

(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]/.θp->θ))[[2]]/(R^(-7/4))/.θ->Sqrt[1.16951]/.m0Sol

h0Sol=Solve[-2 Gls[[3]]==(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(-7/4))/.θ->Sqrt[1.16951]/.m0Sol,h0][[1]]

MCaselle=Table[{hCaselle[θ]((t[θ])^(-Δ[2])),m0 θ/((((θ^2)-1))^(β[2]))}/.m0Sol/.h0Sol,{θ,1.0001,Sqrt[1.16951],0.0001}];

ListPlot[{Mξ3,MCaselle},Joined->True]

χCaselle=Table[{hCaselle[θ]((t[θ])^(-Δ[2])),(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(-7/4))}/.m0Sol/.h0Sol,{θ,1.0001,Sqrt[1.16951],0.0001}];

ListPlot[{χξ3,χCaselle},Joined->True]

χ3Caselle=Table[{hCaselle[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]))[[2]],{{R,θp}}])[[2]](((R t[θ]))^(7/4+Δ[2]))/.θp->θ}/.m0Sol/.h0Sol,{θ,1.0001,Sqrt[1.16951],0.0001}];

ListLogLogPlot[{χ3ξ3,χ3Caselle},Joined->True]