asymptotic expansion of the error function



A workaround might be something similar to (likely this is buggy)

(%i33) erfc_asymp(e, x, n) :=
   subst('erfc = lambda([s], 
     if freeof(x,s) then funmake('erfc, [s]) else
      block([z : taylor(s, x, inf, n)], exp(-s^2) * sum((-1)^k * (2*k+1)!! / (2 * z^(2*k)),k,0,n))), e)$

The stirling function (load("stirling")) works something like this.  Example 

(%i43) erfc_asymp(erfc(x), x, 3);
(%o43)/T/ (1/2-3/(2*x^2)+15/(2*x^4)-105/(2*x^6)+...)*%e^(-x^2)+...

(%i44) erfc_asymp(erfc(x + 1/x), x, 3);
(%o44)/T/ (1/(2*%e^2)-2/(%e^2*x^2)+49/(4*%e^2*x^4)-295/(3*%e^2*x^6)+23161/(48*%e^2*x^8)-101269/(60*%e^2*x^10)+6339529/(1440*%e^2*x^12)-25829911/(2520*%e^2*x^14)+...)*%e^(-x^2)+...

--bw