◼ Input of system  (6.51)

xdot=I  (x-(∑_(k=0)^3)   a[2-k,k] x^(3-k) y^k)

ydot=-I  (y-(∑_(k=0)^3)  b[k,2-k] x^k  y^(3-k))

◼ Normal form (4.22)  of (6.51)  up to the 11th order

x1dot=I   x1 +x1  (∑_(k=1)^5) Y1[2 k ] x1^k y1^k

y1dot=-I   y1+y1 (∑_(k=1)^5) Y2[2  k] x1^k  y1^k

◼ Transformation (4.21)

xsub=x1+(∑_(s=3)^11)(∑_(k=0)^s)(h1[k-1,s-k]   x1^(k)  y1^(s-k) )//Expand;

ysub=y1+(∑_(s=3)^11)(∑_(k=0)^s)(h2[k,s-k-1]  (  x1^(k)  y1^(s-k) ))//Expand;

◼ Choose the distinguishing  transformation

Do[{h1[k,k]=0;h2[k,k]=0},{k,1,6}]

◼ Create equation (2.35)

v1=(D[xsub, x1]   x1dot +D[xsub, y1]   y1dot-(xdot/.{x->xsub,y->ysub})) //Expand;

v2=(D[ysub, x1]   x1dot +D[ysub, y1]   y1dot-(ydot/.{x->xsub,y->ysub})) //Expand;

◼ sh[s]  is a code to solve (2.35) in the space H_s of vector homogeneous functions

sh[s_]:=Module[{g},t1[s-1,0]=First[Solve[Coefficient[v1/.y1->0,x1^s ]==0,h1[s-1,0]]];h1[s-1,0]=h1[s-1,0]/.t1[s-1,0];t1[-1,s]=First[Solve[Coefficient[v1/.x1->0,y1^s ]==0,h1[-1,s]]];h1[-1,s]=h1[-1,s]/.t1[-1,s];Do[If[s-1-i≠i,{t1[s-1-i,i]=First[Solve[Coefficient[v1,x1^(s-i)  y1^i ]==0,h1[s-1-i,i]]];h1[s-1-i,i]=h1[s-1-i,i]/.t1[s-1-i,i]},{t1[s-1-i,i]=First[Solve[Coefficient[v1,x1^(s-i)  y1^i ]==0,Y1[2  i]]];Y1[2 i]=Y1[2  i]/.t1[s-1-i,i]}],{i,1,s-1}];t2[s,-1]=First[Solve[Coefficient[v2/.y1->0,x1^s ]==0,h2[s,-1]]];h2[s,-1]=h2[s,-1]/.t2[s,-1];t2[0,s-1]=First[Solve[Coefficient[v2/.x1->0,y1^s ]==0,h2[0,s-1]]];h2[0,s-1]=h2[0,s-1]/.t2[0,s-1];Do[If[s-1-i≠i,{t2[i,s-1-i]=First[Solve[Coefficient[v2,x1^(i)  y1^(s-i) ]==0,h2[i,s-1-i]]];h2[i,s-1-i]=h2[i,s-1-i]/.t2[i,s-1-i]},{t2[s-1-i,i]=First[Solve[Coefficient[v2,x1^(i)  y1^(s-i) ]==0,Y2[2  i]]];Y2[2 i]=Y2[2  i]/.t2[s-1-i,i]}],{i,1,s-1}]  ]

◼ Calculations using sh

Do[sh[k], {k, 3, 9}]

◼ Output of the coefficients of the normal form

Do[Print[{Y1[2 k],Y2[2 k]}],{k,1,4}]