BeginPackage["Schofield`"]

$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 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
h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] θ^(2i+1), {i, 0, n}]

RFf[y_] := π^-1 y Exp[1/y] ExpIntegralEi[-1/y]
RFnonsingular[n_][θ_] := A RFf[B(θc+θ)] + Sum[FC[i] θ^(2i) , {i, 0, n}]
RF[n_][θ_] := A RFf[B(θc-θ)] + RFnonsingular[n][θ]

RFfCoeff[m_] := 0 /; m < 2
RFfCoeff[m_] := Gamma[m - 1] / π /; m > 1

ruleB[g_] := π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2])

eqLow[n_, h_][m_] := -A B^m RFfCoeff[m] + D[
  RFnonsingular[n][θ] - t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
  {θ, m} ] / m! /. θ -> θc

eqHigh[n_, h_][m_] := D[
  RF[n][θ] - (-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)),
  {θ, m} ] /. θ -> 0

eqMid[n_, h_][m_] := D[
  RF[n][θ] - 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
}

Ghs = {
  0,
  0,
  -1.8452280782328,
  0,
  8.333711750,
  0,
  -95.16896,
  0,
  1457.62,
  0,
  -25891
}

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[{B -> ruleB[g]}, ΦRules, GlRules, GhRules]

eq[F_, g_][m_] := Flatten[Join[{eqLow[F, g][#], eqMid[F, g][#]} & /@ Range[0, m], eqHigh[F, g] /@ Range[0, m, 2]]] /. rules[g]

EndPackage[]