c0=Integrate[u*(1+Erf[u])*Exp[-u^2],{u,-Infinity,Infinity}]
Print[c0]
c1=Integrate[u^2*u*(1+Erf[u])*Exp[-u^2],{u,-Infinity,Infinity}]
Print[c1]

n=20;
mydata=Table[0,{i,1,20}];
For[i=1, i<=n, i++,
  ci=Integrate[u^(i)*u*(1+Erf[u])*Exp[-u^2],{u,-Infinity,Infinity}];
  Print["i=",i,",  ci=",ci];
  mydata[[i]]=ci
]
Print[InputForm[mydata]]
(* Output
{Sqrt[Pi]/2, 5/(4*Sqrt[2]), (3*Sqrt[Pi])/4, 43/(16*Sqrt[2]), (15*Sqrt[Pi])/8, 531/(64*Sqrt[2]), (105*Sqrt[Pi])/16, 8601/(256*Sqrt[2]), 
 
   (945*Sqrt[Pi])/32, 172965/(1024*Sqrt[2]), (10395*Sqrt[Pi])/64, 4161555/(4096*Sqrt[2]), (135135*Sqrt[Pi])/128, 116658675/(16384*Sqrt[2]), 
 
   (2027025*Sqrt[Pi])/256, 3735104625/(65536*Sqrt[2]), (34459425*Sqrt[Pi])/512, 134498225925/(262144*Sqrt[2]), (654729075*Sqrt[Pi])/1024, 
 
   5380583766075/(1048576*Sqrt[2])}

ConditionalExpression[((-1 + (-1)^(2*n))*Gamma[1 + n])/2 - ((1 + (-1)^(2*n))*Gamma[3/2 + n]*Hypergeometric2F1[1/2, 3/2 + n, 3/2, -1])/Sqrt[Pi], 
 Re[n] > -1]


----- for asir
 // DiffValue[I] is int(u^I*sigma(u)*exp(-u^2))
 sqpi=@pi^(1/2); sq2=2^(1/2);
 DiffValue=
 [1/sq2,sqpi/2, 5/(4*sq2), (3*sqpi)/4, 43/(16*sq2), (15*sqpi)/8, 531/(64*sq2), (105*sqpi)/16, 8601/(256*sq2), 
 
   (945*sqpi)/32, 172965/(1024*sq2), (10395*sqpi)/64, 4161555/(4096*sq2), (135135*sqpi)/128, 116658675/(16384*sq2), 
 
   (2027025*sqpi)/256, 3735104625/(65536*sq2), (34459425*sqpi)/512, 134498225925/(262144*sq2), (654729075*sqpi)/1024, 
 
   5380583766075/(1048576*sq2)];

*)

Clear[n]
cn=Integrate[-u^(2*n)*u*(1+Erf[u])*Exp[-u^2],{u,-Infinity,Infinity}]
Print[cn]


