%I #19 Sep 08 2022 08:45:55
%S 0,0,1,2,6,15,46,137,460,1557,5698,21269,83948,339177,1430502,6178981,
%T 27636512,126500209,596320914,2873324677,14203422044,71669915585,
%U 369941778510,1946679921381
%N a(0)=0, a(1)=0; for n>1, a(n) = a(n-1) + (n-1)*a(n-2) + 1.
%H Vincenzo Librandi, <a href="/A186738/b186738.txt">Table of n, a(n) for n = 0..800</a>
%F a(n) = 2*a(n-1)+(n-2)*a(n-2)-(n-2)*a(n-3) with a(0)=a(1)=0, a(2)=1. - _Vincenzo Librandi_, Dec 24 2012
%F E.g.f.: exp(x*(2+x)/2)*sqrt(Pi/2)*(erf(x/sqrt(2))+sqrt(exp(1))*(erf(1/sqrt(2))-erf((1+x)/sqrt(2)))). - _Vaclav Kotesovec_, Dec 26 2012
%F a(n) ~ c * n^(n/2)*exp(sqrt(n)-n/2) * (1+7/(24*sqrt(n))), where c = sqrt(Pi)/2*exp(1/4)*(erf(1/sqrt(2))-1+exp(-1/2)) = 0.32911457143367583724... - _Vaclav Kotesovec_, Dec 26 2012
%t Join[{0}, RecurrenceTable[{a[1] == 0, a[2] == 1, a[3] == 2, a[n] == 2 a[n - 1] + (n - 2) a[n - 2] - (n - 2) a[n - 3]}, a, {n, 40}]] (* _Vincenzo Librandi_, Dec 24 2012 *)
%t FullSimplify[CoefficientList[Series[E^(x*(2+x)/2)*Sqrt[Pi/2]*(Erf[x/Sqrt[2]]+Sqrt[E]*(Erf[1/Sqrt[2]]-Erf[(1+x)/Sqrt[2]])), {x, 0, 20}], x]* Range[0, 20]!] (* _Vaclav Kotesovec_, Dec 26 2012 *)
%o (Magma) I:=[0,0,1,2]; [n le 4 select I[n] else 2*Self(n-1)+(n-3)*Self(n-2)-(n-3)*Self(n-3): n in [1..30]]; // _Vincenzo Librandi_, Dec 24 2012
%K nonn
%O 0,4
%A _Olivier GĂ©rard_, Nov 02 2012
%E Edited by _Bruno Berselli_, Dec 24 2012