Let $ \alpha \ge 0 $ and let $ A \in {\mathbb R} $ and $ B \in {\mathbb R} $ and $ u_* \in {\mathbb R} $. Consider a following integral:
\begin{equation} {\mathfrak I}_{u_*}^{\alpha,A,B}:= \int\limits_{u_*}^\infty u^\alpha \operatorname{erf}(A+B u) \frac{e^{-u^2/2}}{\sqrt{2\pi}} du \end{equation}
We have compute this quantity for in the special case when $\alpha=0$. We have:
\begin{eqnarray} &&{\mathfrak I}_{u_*}^{0,A,B}=\left(\right. \\ && \left. 2 \pi T\left(\frac{A}{\sqrt{B^2+\frac{1}{2}}},\frac{2 A B+2 B^2 u_*+u_*}{\sqrt{2} A}\right)+2 \pi T\left(u_*,\frac{\sqrt{2} (A+B u_*)}{u_*}\right) \right. \\ % && \left. -sign(A \cdot u_*) \cdot \frac{\pi}{2} \right. \\ % && \left. \right)\frac{1}{\pi}+ \frac{1}{2} \text{erf}\left(\frac{A}{\sqrt{2 B^2+1}}\right) \end{eqnarray}
Here $T(,)$ is the Owen's T function.
The result above has been obtained by differentiating the definition with respect to $A$, then doing the integration over $u$ and then integrating the result with respect to $A$ from $A$ to plus infinity and using the result from here. As always, the Mathematica code snippet verifies the result.
In[5169]:=
ll = {};
For[II = 1, II <= 100, II++,
{A, B, us} = RandomReal[{-10, 10}, 3];
x1 = NIntegrate[Erf[A + B u] phi[u], {u, us, Infinity}];
x2 = 1/2 Erf[A /Sqrt[1 + 2 B^2]] +
1/ \[Pi] (-Sign[A us] Pi/2 +
2 \[Pi] OwenT[A/Sqrt[1/2 + B^2], (2 A B + us + 2 B^2 us)/(
Sqrt[2] A)] + 2 \[Pi] OwenT[us, (Sqrt[2] (A + B us))/us]);
ll = Join[ll, {x1 - x2}]
];
ll // Chop
Out[5171]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
5.04842*10^-10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.27066*10^-8, \
0, 0, 0, -1.10895*10^-9, 0, -2.85414*10^-10, 0, 0, 0, 0, 0, 0, 0, 0, \
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
0, 2.08698*10^-10, 0, -1.80319*10^-9, 0, 0, 0, -1.23154*10^-8, 0, 0, \
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
0, 0, 0}
Now, the question is pretty obvious, how do we derive the result in the generic case $\alpha \neq 0 $.
I think I have an answer, although, it isn't pretty...
$${1 \over {\sqrt{2\pi}}}\int_{u_{*}}^{\infty}{u^{\alpha}}{\left(erf(A+Bu)\right)}e^{{-u^2}\over{2}}du$$
$$={\sqrt{2}\over \pi}\int_{u_{*}}^{\infty}{u^{\alpha}}\left(\int_{0}^{A+Bu}e^{-t^2}dt\right)e^{{-u^2}\over{2}}du$$
$$={\sqrt{2}\over \pi}\int_{u_{*}}^{\infty}{u^{\alpha}}\left(\sum_{n=0}^{\infty}{(-1)^n \over n!}\int_{0}^{A+Bu}{t^{2n}}dt\right)e^{{-u^2}\over{2}}du$$
$$={\sqrt{2}\over \pi}\int_{u_{*}}^{\infty}{u^{\alpha}}\left(\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}{(A+Bu})^{2n+1}\right)e^{{-u^2}\over{2}}du$$
So now we have;
$${\sqrt{2}\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\int_{u_{*}}^{\infty}{u^{\alpha}}\cdot {(A+Bu})^{2n+1}\cdot e^{{-u^2}\over{2}}du$$
We can expand ${(A+Bu})^{2n+1}$ thusly;
$${\sqrt{2}\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\int_{u_{*}}^{\infty}{u^{\alpha}}\cdot\sum_{k=0}^{2n+1}{2n+1 \choose k} (A)^{2n+1-k}(Bu)^{k}\cdot e^{{-u^2}\over{2}}du$$
$$={\sqrt{2}\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\sum_{k=0}^{2n+1}{2n+1 \choose k} (A)^{2n+1-k}(B)^{k}\int_{u_{*}}^{\infty}(u)^{k+\alpha}\cdot e^{{-u^2}\over{2}}du$$
Applying the substitution $j={u^2 \over 2}$ gives us;
$$={1\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\sum_{k=0}^{2n+1}{2n+1 \choose k} (A)^{2n+1-k}(B)^{k}\cdot {(\sqrt{2}})^{k+\alpha}\int_{u_{*}^{2}\over2}^{\infty}(j)^{{k+\alpha-1}\over 2}\cdot e^{-j}dj$$
The remaining integral is simply a special case of the upper incomplete gamma function $\small \Gamma(s,x)= \int_{x}^{\infty}{t^{(s-1)}e^{-t}dt}$;
So, finally, we have the following;
$$\mathfrak{I}_{u_{*}}^{{\alpha,A,B}}={A\over \pi}\sum_{n=0}^{\infty}\sum_{k=0}^{2n+1}{(-1)^{n}\cdot {(\sqrt{2}})^{k+\alpha} \cdot (2n)! \over n! \cdot k! \cdot (2n+1-k)!} (A)^{2n-k}(B)^{k}\cdot \Gamma\left({k+\alpha+1\over2},{u_{*}^{2}\over2}\right)$$
I'm not sure if you specifically wanted a final result that made use of Owen's $T$ function but maybe this disaster of a derivation will prove somewhat useful to you...