(PARI) {a(n)=local(A=1+x, B); for(i=1, n, B=subst(A, x, x/(1-x+x*O(x^n)))/(1-x); A=1+3*intformal((B-A)/x)); n!*polcoeff(A, n)}
(PARI) {a(n)=if(n<0, 0, if(n==0, 1, (n-1)!*sum(k=0, n-1, 3*binomial(n, k)*a(k)/k!)))}
(PARI) {a(n)=n!^2*polcoeff(exp(sum(m=1, n, 3*x^m/(m*m!))+x*O(x^n)), n)}
for(n=0, 30, print1(a(n), ", "))
(SageMath)
@CachedFunction
def a(n): return 1 if (n==0) else factorial(n-1)*sum( 3*binomial(n, j)*a(j)/factorial(j) for j in (0..n-1) )