МАТЕМАТИКА даёт такое приближение того же порядка:
PadeApproximant[Exp[x^2] Erfc[x], {x, 0, {1, 2}}]

Думаю получше будет. За ссылку спасибо.
Замена в приведённой Вами формуле, например, коэффициента при иксе сверху эквивалентна формуле

если я правильно сосчитал. И остальные коэффициенты получены некоторым огрублением Паде.
C интегралом-если не удастся найти готовую формулу, то похоже единственный путь-возводить функцию ошибок в квадрат в явном виде, пытаться в двойном ряде сосчитать внутренний ряд из гамма-функций. А если нужен численный результат-то можно сразу от интеграла брать Паде, не считая его.