[Edit]

The generic coefficients in the recursive formula only depend on the polynomials $P_i$. I'll explain when each one is homogeneous of degree $d_i$. The inhomogeneous case is not nearly as clean, so I'd introduce a new variable with substitution polynomial $1$ to reduce it to this case. To find the coefficient of $f_{n_1,\ldots,n_k}$, write $T=\sum_i n_i$ for the total degree. Then in the recursive formula, the smaller coefficient $f_{m_1,\ldots,m_k}$ can only appear when $T = \sum d_im_i$. To find its multinomial and constant parts, write equations for each subscript in terms of the polynomials $P_i$. If $P_i$ has $j_i$ many terms, you will need $j_i$ indices, say $\tau_{i,1}, \ldots, \tau_{i,j_i}$ which sum to $m_i$. These correspond to expanding $P_i^{m_i}$ via the multinomial theorem, which will give you $\binom{m_i}{\tau_{i,1}, \ldots, \tau_{i,j_i}}$ multiplied by terms $c_{i,\alpha_1, \ldots, \alpha_k}^{\tau_{i,h}}$, where the $c_i$s are the coefficients of and the $\alpha_1, \ldots, \alpha_k$ are the exponents of term $h$ in $P_i$. Yikes!

So now we know what the coefficients in the recursive formula should look like, but we need to sum only over the indices $\tau$ which send us to $f_{n_1,\ldots,n_k}$. Each variable gives us one equation, like $n_1 = \sum_{i,h} \tau_{i, h}\alpha_{1,h}$ where the $\alpha$ terms are the exponents again, and $h$ runs from $1$ to $j_i$ for each $i$. These, along with the defining equations $\sum_h\tau_{i,h}=m_i$ and $T = \sum d_im_i$ give us $2k+1$ linear equations ($1$ for each $n$, one for each $m$, and one for the total degree) with $k+1$ dependencies, so you can reduce it to just $k$ equations, and the sum is taken over all non-negative integer solutions to those.

So, for example, the coefficients of your function $G$ are

$$ g_{n_1,n_2,n_3} = \sum \binom{m_1}{\tau_{1,1}, \tau_{1,2}, \tau_{1,3}}\binom{m_3}{\tau_{3,1}, \tau_{3,2}}2^{\tau_{1,2}}g_{m_1,m_2,m_3} $$

Where the sum runs over all non-negative $\tau, m$ such that:

$$ 2m_1 + m_2 + 2m_3 = n_1 + n_2 + n_3 $$ $$ m_1 = \tau_{1,1} + \tau_{1,2} + \tau_{1,3} $$ $$ m_2 = \tau_{2,2} $$ $$ m_3 = \tau_{3,1} + \tau_{3,2} $$ $$ 2\tau_{1,1} + \tau_{1,2} + \tau_{2,1} + 2\tau_{3,1} + \tau_{3,2} = n_1 $$ $$ \tau_{1,2} + \tau_{1,3} + \tau_{3,2} = n_2 \\ $$ $$ \tau_{1,3} = n_3 $$ Sorry I couldn't get the align environment to work, but you get the picture. You can reduce it to just three equations, one of which is constant, which looks a lot nicer, but that's the generic way to do it.

Instead, consider which smaller terms will contribute to the coefficient of $x^n$. For any $m<n$, $x^m$ will contribute if $x^n$ appears in

Finally, I wrote up some code for your function $F$, but I don't have any place to check the values, so although the reasoning about the recurrence is the same, it may have some errors. The coefficients $f_{k,m,n}$ of $x^ky^mz^n$ are given by

$$ f_{k,m,n} = \sum_{(a,b,c)}\binom{a}{n}\binom{a-n}{m-c-n}f_{a,b,c}; $$ where the sum is over all non-negative integers with $2a+b+2c=k+m+n$, so each one is determined by coefficients with strictly smaller total degree (since the monomials $y^m$ never show up). I split the multinomial into two binomials to make some of the restrictions more obvious: for example, if the degree of $z$ is more than the degree of $y$, then the coefficient is 0, because the second binomial always vanishes. This is clear from the given polynomials: the only monomial with a $z$ also contains $y$.

Again, I'm not sure if I implemented this properly, but I'll give some results near $50$ and you can tell me if they look right:

j=60;
F=numpy.zeros((j,j,j), dtype=object); 
F[1,0,0]=1;
F[2,0,0]=1;
F[1,1,0]=1;
F[0,1,1]=1;

for t in range(3,j):                #t is the total degree of the monomial
    for k in range (0, t+1):        #k is the degree of x
        for m in range (0, t+1-k):  #m is the degree of y
            n=t-k-m;                #n is the degree of z
            if m >=n:
                count=0;
                for c in range(0, m-n+1):
                    for a in range (max(n,m-n-c), (t+1)/2):
                        b=t-2*a - 2*c;
                        if b >=0:
                            count=count+binomial(a,n)*binomial(a-n,m-c-n)*F[a,b,c];
                F[k,m,n]=count;
                
for r in range(54, j):
    F[r,0,0]

So from 54-59 inclusive, I got:

I picked that string because they follow some coefficients which about the same size as the one on $x^{54}$, then decrease, then increase again. Seems unlikely to accidentally agree with your function, and may indicate a similar oscillation.

Instead, consider which smaller terms will contribute to the coefficient of $x^n$. For any $m < n$, $x^m$ will contribute if $x^n$ appears in

Finally, I wrote up some code for your function $F$, and checked it against your paper. The coefficients $f_{k,m,n}$ of $x^ky^mz^n$ are given by

$$ f_{k,m,n} = \sum_{(a,b,c)}\binom{a}{n}\binom{a-n}{m-c-n}f_{a,b,c}; $$ where the sum is over all non-negative integers with $2a+b+2c=k+m+n$, so each one is determined by coefficients with strictly smaller total degree (since the monomials $y^m$ never show up). I split the multinomial into two binomials to make some of the restrictions more obvious: for example, if the degree of $z$ is more than the degree of $y$, then the coefficient is 0, because the second binomial always vanishes. This is clear from the given polynomials: the only monomial with a $z$ also contains $y$.

j=60;
F=numpy.zeros((j,j,j), dtype=object); 
F[1,0,0]=1;
F[2,0,0]=1;
F[1,1,0]=1;
F[0,1,1]=1;

for t in range(3,j):                #t is the total degree of the monomial
    for k in range (0, t+1):        #k is the degree of x
        for m in range (0, t+1-k):  #m is the degree of y
            n=t-k-m;                #n is the degree of z
            if m >=n:
                count=0;
                for c in range(0, m-n+1):
                    for a in range (max(n,m-n-c), t/2+1):    #there was a typo here
                        b=t-2*a - 2*c;
                        if b >=0:
                            count=count+binomial(a,n)*binomial(a-n,m-c-n)*F[a,b,c];
                F[k,m,n]=count;

[Edited to fix the above code.]