%I #20 Jan 11 2018 06:25:11
%S 1,1,3,12,76,671,7697,111843,2008907,43535348,1116843468,33429830141,
%T 1153626512945,45418011807581,2021303380235475,100882231471330500,
%U 5607445909955932548,345003931787057849067,23367696786768525050769
%N G.f. A(x) satisfies: A'(x) = 1 + A(x*exp(x)).
%C Higher derivatives of the e.g.f. involve nested exponential functions, e.g.:
%C A''(x)*exp(-x)/(1+x) = 1 + A( x*exp(x)*exp( x*exp(x) ) ).
%F a(n) = Sum_{k=1..n-1} C(n-1,k)* k^(n-k-1)* a(k) for n>1 with a(1)=1.
%e E.g.f.: A(x) = x + x^2/2! + 3*x^3/3! + 12*x^4/4! + 76*x^5/5! +...
%e Related expansions.
%e A'(x) = 1 + x + 3*x^2/2! + 12*x^3/3! + 76*x^4/4! + 671*x^5/5! +...
%e A(x*exp(x)) = x + 3*x^2/2! + 12*x^3/3! + 76*x^4/4! + 671*x^5/5! +...
%t terms = 19; B[_] = 0; Do[B[x_] = 1 + Integrate[B[x], x] /. x -> x Exp[x] + O[x]^terms // Normal, terms];
%t A[x_] = Integrate[B[x], x];
%t CoefficientList[A[x]/x, x] * Range[terms]! (* _Jean-François Alcover_, Sep 15 2011, updated Jan 11 2018 *)
%o (PARI) {a(n)=local(A=x+x^2);for(i=1,n,A=x+intformal(subst(A,x,x*exp(x+1*O(x^n)))));n!*polcoeff(A,n)}
%o (PARI) {a(n)=if(n<1,0,if(n==1,1,sum(k=1,n-1,binomial(n-1,k)*k^(n-k-1)*a(k))))}
%Y Cf. A003659.
%K nonn
%O 1,3
%A _Paul D. Hanna_, Jul 17 2011