nmax = 23;
A087903[n_, k_] :=
A087903[n, k] = StirlingS2[n-1, k] + Sum[ (k-d-1)*
A087903[n-j-1, k-d]*StirlingS2[j, d], {d, 0, k-1}, {j, 0, n-2}]; a[n_] := Sum[
A087903[n, k], {k, 1, n-1}]; a[1] = 1; Table[a[n], {n, 1, nmax}](*
Jean-François Alcover, Oct 04 2011, after
Philippe Deléham *)
Clear[t, n, k, i, nn, x]; coeff = ConstantArray[1, 23]; mp[m_, e_] := If[e==0, IdentityMatrix@ Length@ m, MatrixPower[m, e]]; nn = Length[coeff]; cc = Range[nn]*0 + 1; Monitor[ Do[Clear[t]; t[n_, 1] := t[n, 1] = cc[[n]];
t[n_, k_] := t[n, k] = If[n >= k,
Sum[t[n - i, k - 1], {i, 1, 2 - 1}] +
Sum[t[n - i, k], {i, 1, 2 - 1}], 0];
A4 = Table[Table[t[n, k], {k, 1, nn}], {n, 1, nn}];
A5 = A4[[1 ;; nn - 1]]; A5 = Prepend[A5, ConstantArray[0, nn]];
cc = Total[
Table[coeff[[n]]*mp[A5, n - 1][[All, 1]], {n, 1,
nn}]]; , {i, 1, nn}], i]; cc