■ The operator (3.71) for  system (3.131)

l1[nu1_,nu2_,nu3_,nu4_,nu5_,nu6_]:= 1 nu1 + 0 nu2  -  nu3 + 2 nu4 +1 nu5 + 0 nu6;
l2[nu1_,nu2_,nu3_,nu4_,nu5_,nu6_]:= 0 nu1 + 1 nu2 + 2 nu3  - nu4 + 0 nu5 + 1 nu6;

■ Set  α=β=1  to compute the focus quantities
■ Set  α=1, β=0 to compute the linearizability quantities I_{kk}
■ Set  α=0, β=1 to compute J_{kk}

α=1;  β=1;

■ Definition of  function (4.53)

v[k1_,k2_,k3_,k4_,k5_,k6_]:=v[k1,k2,k3,k4,k5,k6]=
Module[{us,coef},coef=l1[k1,k2,k3,k4,k5,k6]-l2[k1,k2,k3,k4,k5,k6]; us=0; 
v[0,0,0,0,0,0]=1;
If[k1>0,us=us+(l1[k1-1,k2,k3,k4,k5,k6]+α)*v[k1-1,k2,k3,k4,k5,k6]]; 
If[k2>0,us=us+(l1[k1,k2-1,k3,k4,k5,k6]+α)*v[k1,k2-1,k3,k4,k5,k6]];
If[k3>0,us=us+(l1[k1,k2,k3-1,k4,k5,k6]+α)*v[k1,k2,k3-1,k4,k5,k6]];
If[k4>0,us=us-(l2[k1,k2,k3,k4-1,k5,k6]+β)*v[k1,k2,k3,k4-1,k5,k6]];
If[k5>0,us=us-(l2[k1,k2,k3,k4,k5-1,k6]+β)*v[k1,k2,k3,k4,k5-1,k6]];
If[k6>0,us=us-(l2[k1,k2,k3,k4,k5,k6-1]+β)*v[k1,k2,k3,k4,k5,k6-1]];
If [coef!=0,   us=us/coef]; If [coef==0, gg[k1,k2,k3,k4,k5,k6]=us;  us=0];  us]

■ gmax is the number of the focus or  linearizability quantity to be computed

gmax=3;

■ Computing the quantities q[1], q[2], ... up to the order "gmax"

Do[k= sc; num=k;  q[num]=0;
For [i1=0,i1<=2 k,i1++ ,
For[i2=0,i2<=(2 k-i1),i2++ ,
For[i3=0,i3<=(2 k-i1-i2),i3++ ,
For[i4=0,i4<=(2 k-i1-i2-i3),i4++,
For[i5=0,i5<=(2 k-i1-i2-i3-i4),i5++,
For[i6=0,i6<=(2 k-i1-i2-i3-i4-i5),i6++,
If[(l1[i1,i2,i3,i4,i5,i6]==k) &&(l2[i1,i2,i3,i4,i5,i6]==k), v[i1,i2,i3,i4,i5,i6]; 
q[num]=q[num]+gg[i1,i2,i3,i4,i5,i6]TT[i1,i2,i3,i4,i5,i6]]]]]]]],
{sc,1,gmax}]

■ Definition of monomials of system (3.131)

TT[l1_,l2_,l3_,l4_,l5_,l6_]:=a10^l1 a01^l2 a12^l3  b21^l4   b10^l5  b01^l6

■  Computation of  focus or linearizability  quantities

Do[{g[i] = Factor[q[i]], Print[g[i]]}, {i, 1, gmax}]