BeginPackage["Schofield`"] $Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0} β[D_:2] := Piecewise[ { {1/8, D == 2}, {0.326419, D == 3}, {1/2, D == 4}, {β, True} } ] δ[D_:2] := Piecewise[ { {15, D == 2}, {4.78984, D == 3}, {3, D == 4}, {δ, True} } ] α[D_:2] := Piecewise[ { {0, D == 2}, {0.11008, D == 3}, {0, D == 4}, {α, True} } ] Δ[D_:2] := β[D] δ[D] t[θ_] := ((θ)^2 - 1) (* hBasis = Orthogonalize[x^# & /@ Range[0, 20], Function[{f, g}, Integrate[f g (1 - x^2)^2, {x, -1, 1}]]] h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] hBasis[[2 * i + 2]], {i, 0, n}] /. x -> θ / θc *) h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]θ^(2 * i + 1), {i, 0, n}] RFLow[B_, θc_][θ_] := (2 E^(1/( B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] + E^(1/(B (-\[Theta] + \[Theta]c))) (\[Theta] - \[Theta]c) \ ExpIntegralEi[1/(B \[Theta] - B \[Theta]c)] - E^(1/(B \[Theta] + B \[Theta]c)) (\[Theta] + \[Theta]c) ExpIntegralEi[-(1/( B \[Theta] + B \[Theta]c))]) RFHigh[ξ0_][ξ_] := (1+(ξ/ξ0)^2)^(5/6) - 1 RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2(i+1)), {i, 0, n}] ruleB[g_] := B - π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2]) ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - 0.18930 eqLowRHS[F_][m_] := SeriesCoefficient[F[θ], {θ, θc, m}, Assumptions -> Join[{θ < θc, θ > 1}, $Assumptions]] eqLowLHS[h_][m_] :=D[ t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)), {θ, m} ] / m! /. θ -> θc eqLow[F_, h_][m_] := eqLowRHS[F][m] - eqLowLHS[h][m] eqHighRHS[F_][m_] := D[F[θ], {θ, m} ] /. θ -> 0 eqHighLHS[h_][m_] := D[(-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)), {θ, m} ] /. θ -> 0 eqHigh[F_, h_][m_] := eqHighRHS[F][m] - eqHighLHS[h][m] eqMid[F_, h_][m_] := D[ F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η] /. η -> t[θ] / h[θ]^(1 / Δ[2]), {θ, m} ] /. θ -> 1 Φs = { -1.197733383797993, -0.318810124891, 0.110886196683, 0.01642689465, -2.639978 10^-4, -5.140526 10^-4, 2.08856 10^-4, -4.4819 10^-5 } Gls = { 0, -1.3578383417066, -0.048953289720, 0.038863932, -0.068362119, 0.18388370, -0.6591714, 2.937665, -15.61, 96.76, -679, 5340 } Ghs = { 0, 0, -1.8452280782328, 0, 8.333711750, 0, -95.16896, 0, 1457.62, 0, -25891, 0, 502000 } dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)! ΦRules = MapIndexed[dRule[Φ], Φs]; GlRules = MapIndexed[dRule[Gl], Gls]; GhRules = MapIndexed[dRule[Gh], Ghs]; rules[g_] := Join[ΦRules, GlRules, GhRules] (*ξ0 := 0.18930*) gC[0] := 1 tC[0] := 1 eq[F_, g_][m_] := Flatten[Join[{ruleB[g], ruleθ0[g]},{eqLow[F, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[F, g] /@ Range[2, m, 2]]] //. rules[g] newSol[eqs_, oldSol_, newVars_, δ_:0, opts___] := FindRoot[ eqs, Join[ oldSol /. Rule -> List, Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}] ], MaxIterations -> 500, opts ] EndPackage[]