%I #14 Mar 06 2022 08:30:44
%S 1,1,7,16,100,340,2140,9280,60640,311200,2130400,12400000,89094400,
%T 572694400,4314659200,30085907200,237189548800,1771570816000,
%U 14579806451200,115559342963200,990347730035200,8270586336716800,73633536519040000,644303993752499200,5945918623123379200
%N a(0)=1; a(1)=1; for n>1, a(n) = a(n-1) + 3*n*a(n-2).
%F E.g.f. A(x) satisfies the differential equation 6*A(x) + (3*x + 1)*A'(x) - A''(x) = 0, A(0) = 1, A'(0) = 1.
%F E.g.f.: 1 + sqrt(Pi/6) * (1 + 3*x) * exp((1 + 3*x)^2/6) * (erf((1 + 3*x)/sqrt(6)) - erf(1/sqrt(6))).
%F a(n) ~ erfc(1/sqrt(6)) * sqrt(Pi) * 3^(n/2) * exp(sqrt(n/3) - n/2 + 1/12) * n^((n+1)/2) / 2 * (1 + 55/(72*sqrt(3*n)) + 7561/(31104*n)).
%t RecurrenceTable[{a[0]==1, a[1]==1, a[n]==a[n-1] + 3*n*a[n-2]}, a, {n, 0, 20}]
%t nmax = 20; FullSimplify[CoefficientList[Series[1 + Sqrt[Pi/6] * (1 + 3*x) * E^((1 + 3*x)^2/6) * (Erf[(1 + 3*x)/Sqrt[6]] - Erf[1/Sqrt[6]]), {x, 0, nmax}], x] * Range[0, nmax]!]
%K nonn,easy
%O 0,3
%A _Vaclav Kotesovec_, Feb 13 2022, following a suggestion from _John M. Campbell_