Integral of a shifted Gaussian distribution with an error function In the course of computing a convolution of two functions, I have simplified it to a single variable integral of the form
$$\int_0^\infty xe^{-ax^2+bx}\mathrm{erf}(cx+d) dx$$
where $\mathrm{erf}$ is the error function defined as $$\ \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt.$$
I've looked through A Table of Integrals of the Error Functions, especially the formulas on pages 8 and 9. There are a lot of integrals that are similar, but I couldn't find a way to simplify this integral into something that was in that table. I also have tried differentiating under the integral sign using the constants $c$ and $d$ as parameters, but that only seemed to complicate the integration on the parameter after computing the integral with respect to $x$.
Is there any way to find a closed form of this integral?
 A: The result reads:
\begin{eqnarray}
&&\int\limits_0^\infty x e^{-a x^2-b x} \text{erf}(c x+d) dx= \frac{e^{\frac{b^2}{4 a}}}{4 a} \left( \right.\\
&&  \frac{2 b \left(\arctan\left(\frac{\sqrt{a} (b+2 c d)}{2 a d-b c}\right)+\arctan\left(\frac{2 \sqrt{a} d}{b}\right)-\arctan\left(\frac{c}{\sqrt{a}}\right)\right)}{\sqrt{\pi } \sqrt{a}} +\\
&& 2 e^{-\frac{b^2}{4 a}} \text{erf}(d)-\frac{\sqrt{\pi } b \text{erf}\left(\frac{2 a d-b c}{2 \sqrt{a \left(a+c^2\right)}}\right)}{\sqrt{a}}+\frac{2 c e^{-\frac{(b c-2 a d)^2}{4 a \left(a+c^2\right)}}
   \text{erfc}\left(\frac{b+2 c d}{2 \sqrt{a+c^2}}\right)}{\sqrt{a+c^2}} +\\
&&  -\frac{4 \sqrt{\pi } b \left(T\left(\frac{2 a d-b c}{\sqrt{2} \sqrt{a \left(a+c^2\right)}},\frac{\sqrt{a} (b+2 c d)}{2 a d-b c}\right)+T\left(\frac{b}{\sqrt{2} \sqrt{a}},\frac{2 \sqrt{a}
   d}{b}\right)\right)}{\sqrt{a}}\\
&&\left.\right)
\end{eqnarray}
where $T(\cdot,\cdot)$ is the Owen's t function https://en.wikipedia.org/wiki/Owen%27s_T_function .
For[count = 1, count <= 200, count++,
  {a, b, c} = RandomReal[{0, 5}, 3, WorkingPrecision -> 50];
  d = RandomReal[{-5, 5}, WorkingPrecision -> 50];
  I1 = NIntegrate[x Exp[-a x^2 - b x] Erf[c x + d], {x, 0, Infinity}, 
    WorkingPrecision -> 20];
  1/a NIntegrate[
    x Exp[-x^2 - b/Sqrt[a] x] Erf[c/Sqrt[a] x + d], {x, 0, Infinity}];
  1/a Exp[ 
    b^2/(4 a)] NIntegrate[(x - b/(2 Sqrt[a])) Exp[-x^2 ] Erf[
      c/Sqrt[a] ( x - b/(2 Sqrt[a])) + d], {x, b/(2 Sqrt[a]), 
     Infinity}];

  1/a Exp[ 
    b^2/(4 a)] (NIntegrate[(x) Exp[-x^2 ] Erf[
        c/Sqrt[a] ( x - b/(2 Sqrt[a])) + d], {x, b/(2 Sqrt[a]), 
       Infinity}] - 
     b Sqrt[Pi]/(2 Sqrt[a]) 2 T[b/Sqrt[2 a], 
       c/Sqrt[a], (-((b c)/(2 a)) + d) Sqrt[2]]);
  I2 = E^(b^2/(4 a))/(
    4 a) ((2 b)/(
       Sqrt[a] Sqrt[\[Pi]]) (-ArcTan[c/Sqrt[a]] + 
         ArcTan[(2 Sqrt[a] d)/b] + 
         ArcTan[(Sqrt[a] (b + 2 c d))/(-b c + 2 a d)]) + 
      2 E^(-(b^2/(4 a))) Erf[d] - (b Sqrt[\[Pi]])/Sqrt[a]
        Erf[(-b c + 2 a d)/(2 Sqrt[a (a + c^2)])] + (2 c)/Sqrt[
       a + c^2] E^(-((b c - 2 a d)^2/(4 a (a + c^2))))
        Erfc[(b + 2 c d)/(2 Sqrt[a + c^2])] - (4 b Sqrt[\[Pi]])/Sqrt[
       a] (OwenT[b/(Sqrt[2] Sqrt[a]), (2 Sqrt[a] d)/b] + 
         OwenT[(-b c + 2 a d)/(Sqrt[2] Sqrt[a (a + c^2)]), (
          Sqrt[a] (b + 2 c d))/(-b c + 2 a d)]));
  If[Abs[I2/I1 - 1] > 10^(-3), 
   Print["results do not match", {a1, a2, c, d, {I1, I2}}]; Break[]];
  If[Mod[count, 10] == 0, PrintTemporary[count]];
  ];

A: Some remarks. Too long for a comment.
Since our integral can be obtained by differentiation wrt $b$ of a slightly more simple integral, let's consider it:
$$g(a,b,c,d)=\int_0^\infty e^{-ax^2+bx}\mathrm{erf}(cx+d) dx$$
First, some scaling relations:
$$g(a,b,c,d)=\frac{1}{c} g \left(\frac{a}{c^2},\frac{b}{c},1,d \right)=\frac{1}{b} g \left(\frac{a}{b^2},1,\frac{c}{b},d \right)=\frac{1}{\sqrt{a}} g \left(1,\frac{b}{\sqrt{a}},\frac{c}{\sqrt{a}},d \right)$$
Second, by integrating by parts we have an amuzing functional equation:
$$\sqrt{a} e^{-\frac{b^2}{4a}} g(a,b,c,d)+c e^{-d^2} g \left(c^2,-2cd,\sqrt{a},  \frac{b}{2\sqrt{a}}\right)=\frac{\sqrt{\pi}}{2} \left(1-\mathrm{erf} \left(-\frac{b}{2\sqrt{a}} \right) \mathrm{erf}(d) \right)$$
Not sure if this is of any use, but maybe it could suggest values for some particular cases.
By the way, Mathematica refuses to provide explicit form even for $g(1,1,1,1)$ so I wouldn't expect anything for the more general case.
A: Not a full answer, more a long comment, but it might get you started. Define $$f_n:=\int_0^\infty x^n e^{-ax^2+bx}erf(cx+d)dx$$so we want $f_1=-\frac{1}{b}\partial_b f_0$. Defining $u:=\text{erf}(cx+d),\,v:=e^{-ax^2+bx}$, integration by parts gives $$(b+\frac{2a}{b}\partial_b)f_0=-2af_1+bf_0=-\text{erf}(d)-\frac{2c}{\sqrt{\pi}}\int_0^\infty e^{-ax^2+bx-(cx+d)^2}dx.$$The right-hand side is easy enough to write in terms of the error function. 
