with(combinat);

pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;

    res := ind;

    polyvars := indets(poly);
    indvars := indets(ind);

    for v in indvars do
        pot := op(1, v);

        subs1 :=
        [seq(polyvars[k]=polyvars[k]^pot,
             k=1..nops(polyvars))];

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;

pet_cycleind_symm :=
proc(n)
local p, s;
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

V :=
proc(n, k)
    option remember;
    local base, comb, gf, var;

    base := add(Q[p], p=1..k);

    comb := pet_varinto_cind(base, pet_cycleind_symm(n));
    gf :=
    expand(pet_varinto_cind(comb, pet_cycleind_symm(k)));

    for var in [seq(Q[p], p=1..k)] do
        gf := coeff(gf, var, n);
    od;

    gf;
end;

seq(seq(V(d-n, n), n=0..d), d=0..10);