There is another approximation for erf(x) function: erf(x)^2=1-exp(-x^2*(4/Pi+a*x^2)/(1+a*x^2)). Maybe it can be more usefull? At least it's not in series.

https://wikimedia.org/api/rest_v1/media/math/render/svg/caf4121b6bc043fe1dd0138a2e2ade49ee90e4d8

Thanks Andrey! I will try it...

Dear Hugo Hernandez,

I agree with Dr. Osetrov and just want to note you, that erfс(x) is very close to Heaviside step-funtion, so your integral is close to \int exp(x^2) dx 

Hi Dr. Hernandez,

I believe that the analytical solution is given by

∫ {exp(x2)*erfc(x)} dx = [π*erfi(x) − 2*x2 2F2(1, 1; 3/2, 2; x2)] / [2*√(π)] + constant

where pFq is the Generalized Hypergeometric Function.

Dear Yew-Chung Chak,

Please, specify, how you evaluate it? Do you use computer algebra system?

Look at the attach, please

Hi Daniil,

I'm unsure whether or not this is the correct approach, but this is how I evaluated the integral after recognizing that the integral involving exponential and error functions has no elementary antiderivatives:

∫ {exp(x2)*erfc(x)} dx 

• = ∫ {exp(x2)*[1 − erf(x)]} dx
• = ∫ {exp(x2) − exp(x2)*erf(x)} dx
• = ∫ {exp(x2) dx − ∫ {exp(x2)*erf(x)} dx
• = (√(π)/2)*erfi(x) − ∫ {exp(x2)*[(2/√(π)) exp(−x2) Σ_{k=0}^{∞} (2k/(2k + 1)!!) x2k+1]} dx
• = (√(π)/2)*erfi(x) − (2/√(π))*∫ [ Σ_{k=0}^{∞} (2k/(2k + 1)!!) x2k+1]} dx
• = (√(π)/2)*erfi(x) − (2/√(π))*[x22F2(1, 1; 3/2, 2; x2) / 2]
• = (√(π)/2)*erfi(x) − [x22F2(1, 1; 3/2, 2; x2) / √(π)]
• = [π*erfi(x) − 2x22F2(1, 1; 3/2, 2; x2)] / [2*√(π)] + constant

where pFq is the Generalized Hypergeometric Function.

I like the procedure of Yew-Chung Chak, but also you can find the solution using wolfram alpha  

∫ {exp(x2)*erfc(x)} dx 

= ∫ {exp(x2)*[1 − erf(x)]} dx

= ∫ {exp(x2) − exp(x2)*erf(x)} dx

= ∫ {exp(x2) dx − ∫ {exp(x2)*erf(x)} dx

=(Sqrt[Pi] Erfi[x])/2 - (x^2 _2F2(1, 1; 3/2, 2; x2)/Sqrt[Pi]

Thank you Yew-Chung and Ruben Dario for the hypergeometric solution. I think this is the solution that I was looking for.

Generally with combinations involving the error function one finds a tabulation of results from polynomials, which technically deviate from an exact answer only by another error function. The error decreases as the number of terms increase, and can be made as small as required by the situation. It brings to mind the definitions of differentials and integrals where a limit is reached beyond which further precision adds nothing substantial. You can fit your function to a polynomial for integration then compare the results to series expansions of symbolic functions.

There is no reason why you couldn't create a new symbolic function to be the answer, but it creates an obligation to defend the result. 

-\frac{x^2 \, _2F_2\left(1,1;\frac{3}{2},2;x^2\right)}{\sqrt{\pi }}-\frac{x \Gamma

\left(\frac{1}{2},-x^2\right)}{2 \sqrt{-x^2}}