I understand from the OP that the motivation for this question is to find a series expansion in powers of $t$ of 
$$I(t)=\int _{0}^{t}\exp(-x^2 \operatorname{erf}(x))dx.$$
This follows directly from the series expansion of the error function,
$$I(t)=t-\frac{t^4}{2 \sqrt{\pi }}+\frac{t^6}{9 \sqrt{\pi }}+\frac{2 t^7}{7 \pi }-\frac{t^8}{40 \sqrt{\pi }}-\frac{4 t^9}{27 \pi }+\frac{(\pi -28) t^{10}}{210 \pi ^{3/2}}+O\left(t^{11}\right).$$
Here is a comparison of the full integral (blue) with the series expansion to 10-th order (gold).
<IMG SRC="https://ilorentz.org/beenakker/MO/errorfunction.png"/>