(* generic hamiltonian *)
m1=1;
m2=1;
l1=1;
l2=1;
g=1;



h[x1_,x2_,p1_,p2_]:=((1/(2*l1^2*l2^2*(m2^2*(Sin[x2])^2+m1*m2)))*
                    (p1^2*(m2*l2^2)-
                     2*p1*p2*(m2*l2^2+m2*l1*l2*Cos[x2])+
                     p2^2*((m1+m2)*l1^2+m2*l2^2+2*m2*l1*l2*Cos[x2]))
                  - g*((m1+m2)*l1*Cos[x1]+m2*l2*Cos[x1+x2])) ;
                  h[x1,x2,p1,p2]

v1[t_]:=D[h[x1,x2,p1,p2],p1]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};
v2[t_]:=D[h[x1,x2,p1,p2],p2]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};
w1[t_]:=-D[h[x1,x2,p1,p2],x1]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};
w2[t_]:=-D[h[x1,x2,p1,p2],x2]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};

acgoal = 7; precgoal = 12; x1in = 1.5; x2in = 0.00; p1ini = 0. ; 
    p2ini = 0. ;
tmax = 50;
aa = 
        NDSolve[{x1'[t] == v1[t], p1'[t] == 
              w1[t], x2'[t] == v2[t], p2'[t] == w2[t], 
            x1[0] == x1in, p1[0] == p1ini, x2[0] == x2in, 
            p2[0] == p2ini}, {x1, p1, x2, p2}, {t, 0, tmax}, 
          MaxSteps -> 40000, AccuracyGoal -> acgoal, 
          PrecisionGoal -> precgoal];
ParametricPlot[{(x1[t] /. aa)[[1]], (p1[t] /. 
                aa)[[1]]}, {t, 0, tmax}, PlotPoints -> 2000, 
        AxesLabel -> {x1, p1}, Frame -> True];

eo = h[x1in, x2in, p1ini, p2ini]

\(Plot[
      h[\((x1[t] /. aa)\)[\([1]\)], \((x2[t] /. aa)\)[\([1]\)], \((p1[t] /. 
                aa)\)[\([1]\)], \((p2[t] /. aa)\)[\([1]\)]] - eo, 
	{t, 0, 
        tmax}]\)

ParametricPlot[{(x1[t]/.aa)[[1]],(x2[t]/.aa)[[1]]},{t,0,tmax},
PlotPoints->2000,AxesLabel->{x1,p1},Frame->True];

h[0, 0, 0, 0]

en=7;	ensec=h[x1,0,p1,p2];Clear[ener];
p2sol=p2/.Solve[ensec==ener,p2][[2]];p2solm=p2sol/.{x1->0};p1max=(10*en+30)^0.5;
llmax=20;

Do[
x1in=0.;x2in=0;p1ini=p1max-(p1max/llmax)*ll;p2ini=p2solm/.{p1->p1ini,ener->en};
tmax=205;acgoal=8;precgoal=15;
aa=NDSolve[{x1'[t]==v1[t],p1'[t]==w1[t],
x2'[t]==v2[t],p2'[t]==w2[t],x1[0]==x1in,p1[0]==p1ini,
x2[0]==x2in,p2[0]==p2ini},{x1,p1,x2,p2},
	{t,0,tmax},MaxSteps->2400000, AccuracyGoal->acgoal,PrecisionGoal->precgoal];
	
(* Poincare sections at x2=0, v2>0 *)
step=0.05; ret={};lm=Floor[tmax/step]-1;
Do[xa=Mod[(x2[l*step]/.aa)[[1]],2*Pi];
	If[xa>Pi,xa=xa-2*Pi];
	xb=Mod[(x2[(l+1)*step]/.aa)[[1]],2*Pi];
	If[xb>Pi,xb=xb-2*Pi];
		pa=(v2[l*step]/.aa)[[1]];		
	temp1=Sign[xa]*Sign[xb];
	(* If the variable changes sign, we interpolate linearly
	to find the time to where the condition is satisfied *)
If[temp1<0 && pa>0	, to=step*(l-xa/(xb-xa));
	ret=Union[ret,{to}]],{l,1,lm}];		
di=Dimensions[ret][[1]];
(* we plot x1 and p1 at retun times*)

plo[ll]=ListPlot[taaa=Table[{If[Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]>Pi,
Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]-2*Pi,
Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]],
(p1[ret[[l]]]/.aa)[[1]]},
			{l,1,di-1}],Frame->True,PlotRange->All,PlotStyle->{PointSize[0.015],
			RGBColor[1-ll/(2*llmax),0,ll/(2*llmax)]},
			DisplayFunction->Identity],{ll,1,2*llmax-1}];
			psec=Show[Table[plo[ll],{ll,1,2*llmax-1}],DisplayFunction->$DisplayFunction,
			PlotRange->{{-Pi,Pi},All},
			AspectRatio->1];

x1in = 0. ; x2in = 0; p1ini = -2; 
    p2ini = p2solm /. {p1 -> p1ini, 
          ener -> en};
tmax = 35; acgoal = 8; precgoal = 15;
aa = 
        NDSolve[{x1'[t] == v1[t], p1'[t] == 
              w1[t], x2'[t] == v2[t], p2'[t] == w2[t], 
            x1[0] == x1in, p1[0] == p1ini, x2[0] == x2in, 
            p2[0] == p2ini}, {x1, p1, x2, p2}, {t, 0, tmax}, 
          MaxSteps -> 2400000, AccuracyGoal -> acgoal, 
          PrecisionGoal -> precgoal];

\(\(\(
\)\(ParametricPlot[{\(-Pi\) + 
            Mod[\((x1[t] /. aa)\)[\([1]\)] + Pi, 
              2*Pi], \((p1[t] /. aa)\)[\([1]\)]}, {t, 0, tmax}, 
        PlotPoints -> 2000, AxesLabel -> {x1, p1}, 
        Frame -> True];\)\)\)

ParametricPlot[{-Pi + 
            Mod[(x1[t] /. aa)[[1]] + Pi, 2*Pi], -Pi + 
            Mod[+Pi + (x2[t] /. aa)[[1]], 2*Pi]}, {t, 0, tmax}, 
        PlotPoints -> 2000, AxesLabel -> {x1, p1}, 
        Frame -> True];

Do[
      plo[qq] = 
        Plot[0.03*qq, {x, 0, 1}, 
          PlotStyle -> {RGBColor[1 - 0.1*qq, 0, 0.1*qq], 
              Thickness[0.065]}, DisplayFunction -> Identity], {qq, 1, 
        10}]; Show[Table[plo[qq], {qq, 1, 10}], 
      DisplayFunction -> $DisplayFunction]

en=7;	ensec=h[x1,0,p1,p2];Clear[ener];
p2sol=p2/.Solve[ensec==ener,p2][[2]];p2solm=p2sol/.{x1->0};p1max=(10*en+30)^0.5;
llmax=5;

Do[
x1in=0.;x2in=0;p1ini=-2-0.1*ll;p2ini=p2solm/.{p1->p1ini,ener->en};
tmax=205;acgoal=8;precgoal=15;
aa=NDSolve[{x1'[t]==v1[t],p1'[t]==w1[t],
x2'[t]==v2[t],p2'[t]==w2[t],x1[0]==x1in,p1[0]==p1ini,
x2[0]==x2in,p2[0]==p2ini},{x1,p1,x2,p2},
	{t,0,tmax},MaxSteps->2400000, AccuracyGoal->acgoal,PrecisionGoal->precgoal];
	
(* Poincare sections at x2=0, v2>0 *)
step=0.05; ret={};lm=Floor[tmax/step]-1;
Do[xa=Mod[(x2[l*step]/.aa)[[1]],2*Pi];
	If[xa>Pi,xa=xa-2*Pi];
	xb=Mod[(x2[(l+1)*step]/.aa)[[1]],2*Pi];
	If[xb>Pi,xb=xb-2*Pi];
		pa=(v2[l*step]/.aa)[[1]];		
	temp1=Sign[xa]*Sign[xb];
	(* If the variable changes sign, we interpolate linearly
	to find the time to where the condition is satisfied *)
If[temp1<0 && pa>0	, to=step*(l-xa/(xb-xa));
	ret=Union[ret,{to}]],{l,1,lm}];		
di=Dimensions[ret][[1]];
(* we plot x1 and p1 at retun times*)

plo[ll]=ListPlot[taaa=Table[{If[Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]>Pi,
Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]-2*Pi,
Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]],
(p1[ret[[l]]]/.aa)[[1]]},
			{l,1,di-1}],Frame->True,PlotRange->All,PlotStyle->{PointSize[0.015],
			RGBColor[1-ll/llmax,0,ll/llmax]},
			DisplayFunction->Identity],{ll,0,llmax}];
			psec=Show[Table[plo[ll],{ll,0,llmax}],DisplayFunction->$DisplayFunction,
			PlotRange->{{-Pi,Pi},{-3,-2}},
			AspectRatio->1];

en=7;	ensec=h[x1,0,p1,p2];Clear[ener];
p2sol=p2/.Solve[ensec==ener,p2][[2]];p2solm=p2sol/.{x1->0};p1max=(10*en+30)^0.5;
llmax=5;

Do[
x1in=0.;x2in=0;p1ini=-2.35-0.005*ll;p2ini=p2solm/.{p1->p1ini,ener->en};
tmax=205;acgoal=8;precgoal=15;
aa=NDSolve[{x1'[t]==v1[t],p1'[t]==w1[t],
x2'[t]==v2[t],p2'[t]==w2[t],x1[0]==x1in,p1[0]==p1ini,
x2[0]==x2in,p2[0]==p2ini},{x1,p1,x2,p2},
	{t,0,tmax},MaxSteps->2400000, AccuracyGoal->acgoal,PrecisionGoal->precgoal];
	
(* Poincare sections at x2=0, v2>0 *)
step=0.05; ret={};lm=Floor[tmax/step]-1;
Do[xa=Mod[(x2[l*step]/.aa)[[1]],2*Pi];
	If[xa>Pi,xa=xa-2*Pi];
	xb=Mod[(x2[(l+1)*step]/.aa)[[1]],2*Pi];
	If[xb>Pi,xb=xb-2*Pi];
		pa=(v2[l*step]/.aa)[[1]];		
	temp1=Sign[xa]*Sign[xb];
	(* If the variable changes sign, we interpolate linearly
	to find the time to where the condition is satisfied *)
If[temp1<0 && pa>0	, to=step*(l-xa/(xb-xa));
	ret=Union[ret,{to}]],{l,1,lm}];		
di=Dimensions[ret][[1]];
(* we plot x1 and p1 at retun times*)

plo[ll]=ListPlot[taaa=Table[{If[Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]>Pi,
Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]-2*Pi,
Mod[(x1[ret[[l]]]/.aa)[[1]],2*Pi]],
(p1[ret[[l]]]/.aa)[[1]]},
			{l,1,di-1}],Frame->True,PlotRange->All,PlotStyle->{PointSize[0.015],
			RGBColor[1-ll/llmax,0,ll/llmax]},
			DisplayFunction->Identity],{ll,0,llmax}];
			psec=Show[Table[plo[ll],{ll,0,llmax}],DisplayFunction->$DisplayFunction,
			PlotRange->{{-0.25,0.25},{-2.4,-2.35}},
			AspectRatio->1];

x1in = 0. ; x2in = 0; p1ini = -2.37; 
    p2ini = p2solm /. {p1 -> p1ini, 
          ener -> en};
tmax = 35; acgoal = 8; precgoal = 15;
aa = 
        NDSolve[{x1'[t] == v1[t], p1'[t] == 
              w1[t], x2'[t] == v2[t], p2'[t] == w2[t], 
            x1[0] == x1in, p1[0] == p1ini, x2[0] == x2in, 
            p2[0] == p2ini}, {x1, p1, x2, p2}, {t, 0, tmax}, 
          MaxSteps -> 2400000, AccuracyGoal -> acgoal, 
          PrecisionGoal -> precgoal];

ParametricPlot[{-Pi + 
            Mod[(x1[t] /. aa)[[1]] + Pi, 2*Pi], -Pi + 
            Mod[+Pi + (x2[t] /. aa)[[1]], 2*Pi]}, {t, 0, tmax}, 
        PlotPoints -> 2000, AxesLabel -> {x1, x2}, 
        Frame -> True];

\(\(\(
\)\(ParametricPlot[{\(-Pi\) + 
            Mod[\((x1[t] /. aa)\)[\([1]\)] + Pi, 
              2*Pi], \((p1[t] /. aa)\)[\([1]\)]}, {t, 0, tmax}, 
        PlotPoints -> 2000, AxesLabel -> {x1, p1}, 
        Frame -> True];\)\)\)