The  quadratic system to normalize (system (1)):

xdot=x(1-a10 x-a01  y-a12  x^(-1) y^2)//Expand
ydot=-y(1-b01 y-b10 x -b21 y^(-1) x^2)//Expand

define lowest degrees of \mathcal{Q}:

deg={{-1,2},{0,1},{1,0},{2,-1},{3,-1},{4,-2},{5,-2},{6,-3}}

define a subset of \mathcal{Q} ordered lexicographically:

ls={};Do[ Do[{u={deg[[i]][[1]],deg[[i]][[2]]+j}; AppendTo[ls,u]},{j,0,8}],{i,1,Length[deg]}]

ls

reduce the set using a  bound on monomials:

lis={}; Do[{Q=ls[[i]];If[(Q[[1]]<7&& Q[[2]]<7), AppendTo[lis,Q]]; },{i,1,Length[ls]}]

lis

to get a subset qq of \mathcal{Q}  add {0,0} to lis

qq={{-1,2},{-1,3},{-1,4},{-1,5},{-1,6},{0,0},{0,1},{0,2},{0,3},{0,4},{0,5},{0,6},{1,0},{1,1},{1,2},{1,3},{1,4},{1,5},{1,6},{2,-1},{2,0},{2,1},{2,2},{2,3},{2,4},{2,5},{2,6},{3,-1},{3,0},{3,1},{3,2},{3,3},{3,4},{3,5},{3,6},{4,-2},{4,-1},{4,0},{4,1},{4,2},{4,3},{4,4},{4,5},{4,6},{5,-2},{5,-1},{5,0},{5,1},{5,2},{5,3},{5,4},{5,5},{5,6},{6,-3},{6,-2},{6,-1},{6,0},{6,1},{6,2},{6,3},{6,4},{6,5}}

monomials of the system

sl={{-1,2},{0,0},{0,1},{1,0},{2,-1}}

F={};Do[AppendTo[F,{0,0}],{i,1,Length[qq]}];

F

the coefficients of    system of system (1)

sys={{-a12,0},{1,-1},{-a01,b01},{-a10,b10},{0,b21}}

the   lattice vector field (F)_qq to be normalized :

Do[ {i=FirstPosition[qq,sl[[j]]],F[[i]]=sys[[j]]},{j,1,Length[sl]}]

F

qq

m=Length[qq]

comd[uu_,vv_]:={vv[[2]].uu[[1]] vv[[1]]-uu[[2]].vv[[1]] uu[[1]],uu[[2]]+vv[[2]]}

adj[k_,hhV_,ff_]:=Module[{res,rr,j,u,s},{ res={};Do[AppendTo[res,{0,0}],{s,1,Length[ff]}];u=First[FirstPosition[qq,hhV[[2]]]]; Do[If[u+i<k+1,{rr=comd[hhV,{ff[[i]],qq[[i]]}];If[MemberQ[qq,rr[[2]]],{j=FirstPosition[qq,rr[[2]]];res[[j]]=rr[[1]];}]}],{i,1,Length[qq]}]};res]

ads[k_,i_,hk_,fj_]:=Module[{ff,s,u}, res={};Do[AppendTo[res,{0,0}],{u,1,Length[fj]}]; ff=fj;s=1;  While[s<i,{ff=adj[k,hk,ff];res=res+(1/((s)!)) ff};s++];res]

to compute the normal form till the term (4,4):

k1=FirstPosition[qq,{4,4}][[1]]

lambda={1,-1}

af=F

set k1+1 as the first argument in ads

Do[{Q=qq[[s]];If[Q.lambda ≠0, { h={ 1/(Q.lambda) af[[s]],qq[[s]]};af=af+ ads[k1+1,m,h,af]}]},{s,1,m}]

af

eq1=0; Do[eq1=eq1+af[[i]][[1]]   x^qq[[i]][[1]]   y^qq[[i]][[2]],{i,1,m}]; eq1=x eq1;

eq2=0; Do[eq2=eq2+af[[i]][[2]]   x^qq[[i]][[1]]   y^qq[[i]][[2]],{i,1,m}]; eq2=y eq2;

g1k - the coefficient of  x^{k+1} y^k  in  the first equation of the normal form,
g2k - the coefficient of  x^k  y^{k+1}  in  the second equation of the normal form

eq1

g11=Coefficient[eq1,x^2 y]//Expand

g12=Coefficient[eq2,x y^2]//Expand

g21=Coefficient[eq1,x^3 y^2]//Factor

g22=Coefficient[eq2,x^2  y^3]//Factor

g31=Coefficient[eq1,x^4 y^3]//Factor

g32=Coefficient[eq2,x^3  y^4]//Factor