alias( om = RootOf(Z_^2 + Z_ + 1 = 0) );

X1o := 1/2*T*(3*T-2*u+3)*(1616*o*u^2-2394*o*u+268*u^2+1323*T+2646*o-1008*u-1323
)*(1616*o*u^2-2394*o*u+1348*u^2-1323*T+2646*o-1386*u+3969)*(42*T-37*u-21)*(128*
o*u^2+64*u^2+1323*T-1323*o-1323)*(128*o*u^2+64*u^2-1323*T-1323*o)/(3*T+u+3)/(
1616*o*u^2-2394*o*u+268*u^2+1323*T+2646*o+315*u-1323)/(1616*o*u^2-2394*o*u+1348
*u^2-1323*T+2646*o-2709*u+3969)/(42*T-16*u-21)/(128*o*u^2+64*u^2+1323*T-1323*o-\
1323*u-1323)/(128*o*u^2+64*u^2-1323*T-1323*o+1323*u):
Y1o := -1/2*(T-1)*(3*T-2*u-6)*(1616*o*u^2+2394*o*u+1348*u^2+1323*T+2646*o+1386*
u+2646)*(1616*o*u^2+2394*o*u+268*u^2-1323*T+2646*o+1008*u)*(42*T-37*u-21)*(128*
o*u^2+64*u^2+1323*T-1323*o-1323)*(128*o*u^2+64*u^2-1323*T-1323*o)/(3*T+u-6)/(
1616*o*u^2+2394*o*u+1348*u^2+1323*T+2646*o+2709*u+2646)/(1616*o*u^2+2394*o*u+
268*u^2-1323*T+2646*o-315*u)/(42*T-16*u-21)/(128*o*u^2+64*u^2+1323*T-1323*o-\
1323*u-1323)/(128*o*u^2+64*u^2-1323*T-1323*o+1323*u):
T1o := [[-1/9*u^2+89/486*u^3-1, 1/9*u^2+89/486*u^3+2, -29/1701*u^3+3/14*u+1/2],
[(-1913/11907*o-493/3402)*u^3+(-1427/1323*o-730/1323)*u^2+(8/7*o+2/21)*u-2*o+1,
(1913/11907*o+125/7938)*u^3+(-1427/1323*o-697/1323)*u^2+(-8/7*o-22/21)*u-2-2*o,
53/11907*u^3+u^2*(-128/1323*o-64/1323)+2/3*u+o+1], [(1913/11907*o+125/7938)*u^3
+(1427/1323*o+697/1323)*u^2+(-8/7*o-22/21)*u+2*o+3, (-1913/11907*o-493/3402)*u^
3+(1427/1323*o+730/1323)*u^2+(8/7*o+2/21)*u+2*o, 53/11907*u^3+u^2*(128/1323*o+
64/1323)+2/3*u-o], [-1+(14003/7938*o-1361/7938)*u^3+(-o-1/3)*u^2+2/3*u*o, (-\
27193/7938*o+2713/7938)*u^3+(1513/1323*o-256/1323)*u^2+(-22/7*o-12/7)*u-2-2*o,
(6382/3969*o-1195/3969)*u^3+(-1048/1323*o-524/1323)*u^2+(4/3+2/3*o)*u-o], [(-\
103993/2722734*o-368945/2722734)*u^3+(-10784/9261*o-4177/9261)*u^2+(40/21*o+8/
21)*u-2*o+1, (-103993/2722734*o-368945/2722734)*u^3+(10784/9261*o+4177/9261)*u^
2+(40/21*o+8/21)*u+2*o, u^3*(-48446/1361367*o-130561/1361367)+u*(1/2+2/21*o)+1/
2], [(-27193/7938*o+2713/7938)*u^3+(-1513/1323*o+256/1323)*u^2+(-22/7*o-12/7)*u
+2*o+3, 2+(14003/7938*o-1361/7938)*u^3+(o+1/3)*u^2+2/3*u*o, (6382/3969*o-1195/
3969)*u^3+(1048/1323*o+524/1323)*u^2+(4/3+2/3*o)*u+o+1], [(27193/7938*o+14953/
3969)*u^3+(1513/1323*o+1769/1323)*u^2+(22/7*o+10/7)*u-2*o+1, 2+(-14003/7938*o-\
7682/3969)*u^3+(-o-2/3)*u^2+u*(-2/3-2/3*o), (-6382/3969*o-7577/3969)*u^3+(-1048
/1323*o-524/1323)*u^2+(-2/3*o+2/3)*u-o], [-1+(-14003/7938*o-7682/3969)*u^3+(o+2
/3)*u^2+u*(-2/3-2/3*o), (27193/7938*o+14953/3969)*u^3+(-1513/1323*o-1769/1323)*
u^2+(22/7*o+10/7)*u+2*o, (-6382/3969*o-7577/3969)*u^3+(1048/1323*o+524/1323)*u^
2+(-2/3*o+2/3)*u+o+1], [(103993/2722734*o-132476/1361367)*u^3+(10784/9261*o+
6607/9261)*u^2+(-40/21*o-32/21)*u+2*o+3, (103993/2722734*o-132476/1361367)*u^3+
(-10784/9261*o-6607/9261)*u^2+(-40/21*o-32/21)*u-2-2*o, u^3*(48446/1361367*o-\
82115/1361367)+u*(-2/21*o+17/42)+1/2]]:

X1 := subs(o=om,X1o); Y1 := subs(o=om,Y1o); T1 := subs(o=om,T1o);

#
#  Check that [X1,Y1]  is a solution up to $O(u^3)$.
#

for k from 1 to 9 do
  for j from 1 to 3 do
    xx[j] := convert(series( subs( T = T1[k,j], X1 ), u,4),polynom);
    yy[j] := convert(series( subs( T = T1[k,j], Y1 ), u,4),polynom);
  od;
  print( factor(xx[1]-xx[2]), factor(xx[2]-xx[3]), factor(yy[1]-yy[2]), factor(yy[2]-yy[3]) );
od: