c0=Integrate[D[u*(1+Erf[u]),u]*Exp[-u^2],{u,-Infinity,Infinity}]
Print[c0]
c1=Integrate[u^2*D[u*(1+Erf[u]),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)*D[u*(1+Erf[u]),u]*Exp[-u^2],{u,-Infinity,Infinity}];
  Print["i=",i,",  ci=",ci];
  mydata[[i]]=ci
]
Print[InputForm[mydata]]
(* Output
 c0 Sqrt[Pi] は手動で足す.
{3/(2*Sqrt[2]), Sqrt[Pi]/2, 13/(8*Sqrt[2]), (3*Sqrt[Pi])/4, 101/(32*Sqrt[2]), (15*Sqrt[Pi])/8, 1167/(128*Sqrt[2]), (105*Sqrt[Pi])/16,  
   18147/(512*Sqrt[2]), (945*Sqrt[Pi])/32, 356325/(2048*Sqrt[2]), (10395*Sqrt[Pi])/64, 8458245/(8192*Sqrt[2]), (135135*Sqrt[Pi])/128,  
 235344375/(32768*Sqrt[2]), (2027025*Sqrt[Pi])/256, 7504668675/(131072*Sqrt[2]), (34459425*Sqrt[Pi])/512,  
   269651180925/(524288*Sqrt[2]), (654729075*Sqrt[Pi])/1024}

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

----- for asir
 // DiffValue[I] is int(u^I*sigma{1}(u)*exp(-u^2))
 sqpi=@pi^(1/2); sq2=2^(1/2);
 DiffValue=
 [
 3/(2*sq2), sqpi/2, 13/(8*sq2), (3*sqpi)/4, 101/(32*sq2), (15*sqpi)/8, 1167/(128*sq2), (105*sqpi)/16,  
   18147/(512*sq2), (945*sqpi)/32, 356325/(2048*sq2), (10395*sqpi)/64, 8458245/(8192*sq2), (135135*sqpi)/128,  
 235344375/(32768*sq2), (2027025*sqpi)/256, 7504668675/(131072*sq2), (34459425*sqpi)/512,  
   269651180925/(524288*sq2), (654729075*sqpi)/1024
];

*)

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


