%I #13 May 12 2019 18:16:34

%S 1,1,4,8,18,24,52,60,106,135,213,225,397,411,599,719,1001,1019,1533,

%T 1553,2192,2464,3151,3175,4502,4641,5888,6404,8145,8175,11040,11072,

%U 13863,14811,17886,18390,23723,23761,28440,30140,36650,36692,45952,45996,55095,58535,68084,68132,83720,84193

%N a(1) = 1; a(n+1) = Sum_{d|n} sigma(n/d)*a(d), where sigma = sum of divisors (A000203).

%F G.f. A(x) satisfies: A(x) = x * (1 + Sum_{k>=1} sigma(k)*A(x^k)).

%F G.f.: A(x) = Sum_{n>=1} a(n)*x^n = x * (1 + Sum_{i>=1} Sum_{j>=1} sigma(i)*a(j)*x^(i*j)).

%t a[n_] := a[n] = Sum[DivisorSigma[1, (n - 1)/d] a[d] , {d, Divisors[n - 1]}]; a[1] = 1; Table[a[n], {n, 1, 50}]

%t terms = 50; A[_] = 0; Do[A[x_] = x (1 + Sum[DivisorSigma[1, k] A[x^k], {k, 1, terms}]) + O[x]^(terms + 1) // Normal, terms + 1]; Rest[CoefficientList[A[x], x]]

%t a[n_] := a[n] = SeriesCoefficient[x (1 + Sum[Sum[DivisorSigma[1, i] a[j] x^(i j), {j, 1, n - 1}], {i, 1, n - 1}]), {x, 0, n}]; Table[a[n], {n, 1, 50}]

%o (PARI) lista(nn) = { my(va=vector(nn)); va[1] = 1; for (n=2, nn, va[n] = sumdiv(n-1, d, sigma((n-1)/d)*va[d])); va;} \\ _Michel Marcus_, Apr 30 2019

%Y Cf. A000203, A007557, A307794.

%K nonn

%O 1,3

%A _Ilya Gutkovskiy_, Apr 30 2019