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=\sum_{p=1}^\infty c_p t^p.$$
The coefficients $c_p=p^{-1}d_{p-1}$ follow from the series expansion
$e^{-x^2\,{\rm erf}\,x}=\sum_{p=0}^\infty d_p x^p$, resulting in
$$I(t)=\sum_{p=1}^\infty c_p t^p=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).$$
The series $I(1)=\sum_{p=1}^\infty c_p$ seems to converge:

<IMG SRC="https://ilorentz.org/beenakker/MO/errorfunction_3.png"/>

Plot of $I_N=\sum_{p=0}^{N} c_p $ as a function of $N$ up to $N=25$. The value of $I_{25}=0.8162$ agrees with $I(1)=0.816377$ to three decimal places. For $N=50$ the agreement is up to six decimal places.