How to approximate the difference between two functions by an asymptotic expansion I need to approximate this function
$$f(x) = \operatorname{erfc}(x-a)-\operatorname{erfc}(x),$$
where $x,a \in \mathbb R^+, a \approx 12$, for large values of $x$.
For large $x$, $f(x)$ tends to $1-1 = 0$, and when I calculate it, values of $x > 17$ or so give numerical instability (and in some other systems, even at lower values of $x$).
For the calculations I need to do, $f(x)$ is at the denominator of another function whose numerator also gets close to $0$ when $x$ is large (but is numerically stable):
$$g(y) = \frac {e^{-(x-y)^2}} {\operatorname{erfc}(x-a)-\operatorname{erfc}(x)}$$
So I think I need a better estimate of $f(x)$.
For a similar problem, I was advised to use an asymptotic expansion of $\operatorname{erfc}$ for large values of $x$ (like here), and that worked very well indeed.
So I tried it in this case, and while it also seemed to work, it made me wonder whether I am taking the correct approach.
In the other problem, the overall function to approximate was:
$$h(x) = \frac {e^{-x^2}} {\operatorname{erfc}(x)}$$
There weren't many possibilities: just replace $\operatorname{erfc}$ with its asymptotic series (which in that case also cancelled out the exponential at the numerator).
So here I simply plugged in the expression of $\operatorname{erfc}$ as a power series, truncated at $n=4$, both when the argument of $\operatorname{erfc}$ was $x$ and when it was $x-a$. So I wrote the difference of the $\operatorname{erfc}$'s as a difference between two power series.

Is this a legitimate approach, or should I instead look for the asymptotic series expansion of the difference itself? Or even of the reciprocal of the difference, as that's what goes into the above mentioned $g(y)$?

Perhaps the asymptotic series expansion of the difference would converge more rapidly, allowing to use fewer terms?
Any ideas/advice?
Thanks!

EDIT: following up from Gary's suggestion
This is what I find (using $a = 12, n = 4$):

This puzzles me in several different ways.
First, it looks like the $\operatorname{erfc}$ difference (original function, blue curve) is quite stable at large values of $x$, unlike I had previously found (or assumed?). It's not clearly visible from this plot, but the blue curve continues all the way down to the right. So there seems to be no need for a polynomial approximation.
In fact it's on the 'other side', when $x$ becomes negative, and not even very large in absolute value (i.e. where I thought I was 'safe'), that this particular software has trouble with it (you can see it stops being displayed in log scale at about $x \le -5$, meaning it falls to exactly 0 there). R is not behaving very differently, either:
erfc <- function(x) 2 * pnorm(x * sqrt(2),lower=FALSE)
curve(erfc(x-12)-erfc(x), -20, 50, col = "blue", log = "y")


Gary's function (green curve) and the difference between the asymptotic approximations of $\operatorname{erfc}$ (red curve) seem to work exactly like the $\operatorname{erfc}$ difference on the 'large $x$' side, whereas on the 'negative $x$' side only the red curve seems to work (and it was not even meant to do that! It was designed to work for large positive $x$. I guess this expansion has some property that makes it work for negative $x$, too).
And in fact, using the fact that $\operatorname{erfc}(x) = -\operatorname{erfc}(-x)$:
curve(-erfc(12-x)+erfc(-x), -50, 50, col = "blue", log = "y")


it looks like I can get the whole curve using the original $\operatorname{erfc}$ difference, by defining it piecewise, if I can figure out where to split it (at the moment I would guess $a/2$).
I suppose the polynomial expansion would still be useful if one had the possibility to cancel out some exponential factors. Then one could go potentially to very large values of $x$.
Interesting stuff, but it makes me wonder why there isn't a ready-made definition of $\operatorname{erfc}$ in R, for instance, and why its definition via pnorm performs better on one side than on the other.

EDIT 2 after further work
And in fact I see now that $e^{-x^2}$ does cancel out:
$$g(y) = {{e^ {- \left(x-y\right)^2 }}\over{e^ {- \left(x-a\right)^2 }\,P
 \left(x-a\right)-e^ {- x^2 }\,P\left(x\right)}} = {{e^{-y^2+2\,x\,y+a^2}}\over{e^{2\,a\,x}\,P\left(x-a\right)-e^{a^2}
 \,P\left(x\right)}}$$
 A: You may try the following expansion. By http://dlmf.nist.gov/8.4.E6 and http://dlmf.nist.gov/8.7.E2
\begin{align*}
& \operatorname{erfc}(x - a) - \operatorname{erfc}(x) = \frac{1}{{\sqrt \pi  }}\left( {\Gamma \!\left( {\tfrac{1}{2},(x - a)^2 } \right) - \Gamma \!\left( {\tfrac{1}{2},x^2 } \right)} \right)
\\ &
 = \frac{{\mathrm{e}^{ - (x - a)^2 } }}{{(x - a)\sqrt \pi  }}\sum\limits_{n = 0}^\infty  {( - 1)^n \frac{{(2n - 1)!!}}{{2^n (x - a)^{2n} }}\left( {1 - \mathrm{e}^{ - 2ax + a^2 } \sum\limits_{k = 0}^n {\frac{{(2ax - a^2 )^k }}{{k!}}} } \right)} ,
\end{align*}
whenever $\left| {2ax - a^2 } \right| \le (x - a)^2$. In particular, this works when $x > (\sqrt 2  + 2)a\approx 3.414213562\times a$.
A: Another possibility to try is, for $x \to \infty,$
$$\text{erfc}(x-a) \sim \frac{2}{\sqrt{\pi}}e^{-(x-a)^2} \sum_{m=0}^\infty (2x)^{-(m+1)} H_m(a)$$ where the $H_m(a)$ are the (physicist's, and Mathematica's) Hermite polynomials.  Note by the closed form for $H_m(0)$ this returns the well-known asymptotic formula for erfc.  The subtraction in fact can be written
$$\text{erfc}(x-a) \sim \frac{2}{\sqrt{\pi}}e^{-x^2} \sum_{m=0}^\infty (2x)^{-(m+1)} \Big( e^{a(2x-a)} H_m(a) - H_m(0)  \Big).$$
I haven't done much numerical work with this, but for a~17 and x~20, it seems I needed a good number of terms, like 50, so I don't know if this is useful.  It should be noted that the Hermite polynomials can be calculated stably with Clenshaw recursion, if your goal is to put this in a computer routine limited to double precision.
