(PARI) a(n) = n!/6*sum(k=0, n-3, 1/k!);
(PARI) a(n) = n!*sum(k=0, n, binomial(k, 3)/k!);
(PARI) my(N=30, x='x+O('x^N)); concat([0, 0, 0], Vec(serlaplace(x^3/6*exp(x)/(1-x))))
(PARI) my(N=30, x='x+O('x^N)); concat([0, 0, 0], Vec(sum(k=3, N, k!*x^k/(1-x)^(k+1))/6))