I know that the imaginary error function, $\mathrm{Erfi}(x)=(2/\sqrt{\pi})\int_0^x \exp{t^2} \mathrm{d}t$, has the asymptotic expansion given in the answer to this question: Asymptotic order of $\frac{\mathrm{erfi}(\sqrt{x})}{\exp(x)\sqrt{x}}$ .
My question is about deriving this expansion, which I think is based on integration by parts. For example, if we use integration by parts once we can rewrite the function as $\mathrm{Erfi}(x) = \frac{2}{\sqrt{\pi}} \left[\frac{1}{2t}\mathrm{e}^{t^2}|_0^x+\frac{1}{2}\int_0^x\frac{\mathrm{e}^{t^2}}{t^2} \mathrm{d}t\right]$,
and then keep on doing this procedure for the leftover integral. This indeed produces the terms in the expansion mentioned above. However, what happens to the terms evaluated at $t=0$? For example, in the first term of the RHS above, we get a blow-up when we plug $t=0$.
How do we formally justify ignoring these blowing-up terms?
Note that $$ \operatorname{erfi} (\mathrm{i} \infty) = \frac{2}{\sqrt{\pi}} \int \limits_0^{\mathrm{i} \infty} \mathrm{e}^{t^2} \, \mathrm{d}t = \frac{2 \mathrm{i}}{\sqrt{\pi}} \int \limits_0^\infty \mathrm{e}^{-s^2} \, \mathrm{d}s = \mathrm{i} \, . $$ Therefore we can write $$ \operatorname{erfi}(x) = \mathrm{i} + \frac{2}{\sqrt{\pi}} \int \limits_{\mathrm{i} \infty}^x \mathrm{e}^{t^2} \, \mathrm{d} t \, .$$ Now we can integrate by parts without any problems and since $\mathrm{i}$ is completely negligible compared to $\operatorname{erfi}(x)$ in the limit $x \to \infty$ , we obtain the desired expansion (with $(-1)!! \equiv 1$): $$ \operatorname{erfi}(x) \sim \frac{\mathrm{e}^{x^2}}{\sqrt{\pi} \, x} \sum \limits_{k=0}^\infty \frac{(2k-1)!!}{(2 x^2)^k} \, , \, x \to \infty \, . $$