(PARI) seq(max_n) = {my(a = f = vector(max_n), s, D); f[1]=1;
for(j=1, max_n - 1, f[j+1] = 1/j * sum(k=1, j, sumdiv(k, d, d * f[d]) * f[j-k+1]));
for(n=3, max_n, s=0; forpart(P=n, D=Set(P); if(#D==3, s+=f[P[1]]*f[P[2]]*f[P[3]]; next());
if(#D==1, s+= binomial(f[P[1]]+2, 3); next());
if(P[1] == P[2], s += binomial(f[P[1]]+1, 2) * f[P[3]],
s += binomial(f[P[2]]+1, 2) * f[P[1]]), [1, n], [3, 3]); a[n] = s ); a[3..max_n] }; \\
Washington Bomfim, Dec 22 2020