Haminiz = (px^2+py^2)/2+y*px-x*py -(1-mu)/Sqrt[(x-mu)^2+y^2]-mu/Sqrt[(x+1-mu)^2+y^2];
Sgen = px * (mu+rho*Cos[theta]) + (py - mu) * rho*Sin[theta];
Spx = D[Sgen, px];
Spy = D[Sgen,py];
Srho = D[Sgen,rho];
Stheta = D[Sgen,theta];
solsistvecchimomenti=Simplify[Expand[Solve[{prho== Srho , ptheta == Stheta}, {px,py}]]];
sostcoordinate = {{x -> Spx, y -> Spy}};
sostcompleta = Join[ solsistvecchimomenti[[1]] , sostcoordinate[[1]] ];
HamPolarElioc=Expand[Simplify[Haminiz /. sostcompleta]]

TreeForm[ HamPolarElioc ]
Head[ HamPolarElioc ];
Length[ HamPolarElioc ];
Hamsempl = HamPolarElioc;
Hamsempl -=  HamPolarElioc[[5]] + HamPolarElioc[[6]] + HamPolarElioc[[8]];
Hamsempl += - (1-mu)/rho

Ntronc = 5;
Hamsempl /. { prho -> PX, ptheta -> PY+1, rho -> X+1, theta -> Y - 2*Pi/3};
Series[%,{X,0,Ntronc}];
Series[%,{Y,0,Ntronc}];
Normal[%];
HamsemplPXY = ExpandAll[%];
MonomialList[ ExpandAll[%] , {PX,PY,X,Y}]

SvilPotPrim2 = Normal[ Series[ - mu / Sqrt[1 + u] , {u, 0, Ntronc} ] ];
uXYtmp =(2 * rho * Cos[theta] + rho^2) /. { rho -> X+1, theta -> Y - 2*Pi/3};
uXY = ExpandAll[ Normal[ Series[uXYtmp, {X,0,Ntronc}, {Y,0,Ntronc} ] ] ];
parteComplicataXY = Expand[ Normal[ SvilPotPrim2 /. {u -> uXY} ] ] ;
HamTotTmp = HamsemplPXY + parteComplicataXY;

SelPolOm[expr_,gr_]:=Module[{ris = 0, espan, nummon, j, addendo, gradd},
         espan = MonomialList[ ExpandAll[expr] , {PX,PY,X,Y} ];
	nummon = Length[ espan ];
	For[j=1, j≤nummon, j++,
	      addendo = espan[[j]];
               gradd = 0;
	      gradd += Exponent[addendo, PX ];
	      gradd += Exponent[addendo, PY ];
	      gradd += Exponent[addendo, X ];
	      gradd += Exponent[addendo, Y ];
	      If[ gradd == gr, ris += addendo ];
       ];
          ris]
HamPXY = 0;
For[ grpol = 0, grpol≤ Ntronc, grpol++,
HamPXY += SelPolOm[HamTotTmp, grpol]; ]
HamPXY
HamQuad = SelPolOm[HamTotTmp, 2];
matQ = D[ D[ HamQuad, {{PX, PY,X,Y}} ] , {{PX, PY,X,Y}} ];
matJ={{0,0,-1,0},{0,0,0,-1},{1,0,0,0},{0,1,0,0}};
matJQ = matJ . matQ;
CharacteristicPolynomial[ matJQ , lambda ]