For inverse functions like InverseErf[]
, the usual approach for determining their series coefficients is Lagrangian inversion. (That is, what InverseSeries[]
does.)
In particular, Mathematica supports the (partial) Bell polynomials that appear in the series coefficients for the inverse series as BellY[]
. (See also the discussion in Charalambides's book.) Using this, along with the general formula for the derivative of Erf[]
(formula 7.1.19 in Abramowitz and Stegun), we have
SetAttributes[inverseErfCoefficient, Listable];
inverseErfCoefficient[0] = 0; inverseErfCoefficient[1] = Sqrt[π]/2;
inverseErfCoefficient[m_Integer?Positive] :=
(Sqrt[π]/2)^m BellY[Table[{(m + k - 2)!, (-2)^k/(2 k) Sqrt[π]/Gamma[1 - k/2]},
{k, 2, m}]]/(m! (m - 1)!)
Test:
inverseErfCoefficient[Range[0, 7]]
{0, Sqrt[π]/2, 0, π^(3/2)/24, 0, (7 π^(5/2))/960, 0, (127 π^(7/2))/80640}
Sum[inverseErfCoefficient[k] x^k, {k, 0, 20}] - (InverseErf[x] + O[x]^21)
O[x]^21
As it turns out, the series coefficients for the Taylor expansion of the inverse error function at $0$ satisfy a simple nonlinear recurrence. The formulae used in the following are adapted from Carlitz and Dominici:
ck[0] = 0; ck[1] = 1;
ck[n_Integer?Positive] := ck[n] = 2 Sum[Binomial[n - 1, k] ck[k - 1] ck[n - k], {k, n - 1}]
SetAttributes[inverseErfCoefficient, Listable];
inverseErfCoefficient[m_Integer?NonNegative] := (Sqrt[π]/2)^m ck[m]/m!
This version should give the same results.