Solution (source code)

= Solution

Let
$$
I(x)=\int_x^\infty e^{-t^2}\,dt.
$$
Since $d(e^{-t^2})=-2te^{-t^2}dt$, repeated <integration by parts> gives
$$
\begin{aligned}
I(x)
&=\frac{e^{-x^2}}{2x}
-\frac12\int_x^\infty t^{-2}e^{-t^2}\,dt\\
&=\frac{e^{-x^2}}{2x}
-\frac{e^{-x^2}}{4x^3}
+\frac34\int_x^\infty t^{-4}e^{-t^2}\,dt\\
&=\frac{e^{-x^2}}{2x}
\left(1-\frac1{2x^2}+\frac3{4x^4}
+O(x^{-6})\right).
\end{aligned}
$$
The final remainder estimate follows by one further integration by parts, or by bounding the remaining positive integral by its leading boundary term. Therefore the requested first three terms are
$$
\boxed{\operatorname{erfc}(x)
\sim\frac{e^{-x^2}}{2\sqrt{2\pi}\,x}
\left(1-\frac1{2x^2}+\frac3{4x^4}\right).}
$$