(* Double pendulum. x1in is the initial value for the first angle 
	and x2in for the second angle (see lect. notes) ,
	p1ini and p2ini are related to the initial velocities, tmax is the time 
	you run the ODE *)
x1in=1.5;x2in=0.00;p1ini=0.;p2ini=0.;
tmax=50;


(* you do not need to modify anything below, unless you want to
spend the time to understand it *)
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]};                
vx1[x1_,x2_,p1_,p2_]:=D[h[x1,x2,p1,p2],p1];
vx2[x1_,x2_,p1_,p2_]:=D[h[x1,x2,p1,p2],p2];
vp1[x1_,x2_,p1_,p2_]:=-D[h[x1,x2,p1,p2],x1];
vp2[x1_,x2_,p1_,p2_]:=-D[h[x1,x2,p1,p2],x2]; 
acgoal=7;precgoal=12;
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];

\(\(\( (*Motion\ of\ m2\ *) \)\(
\)\(ParametricPlot[{Sin[\((x1[t] /. 
                    aa)\)[\([1]\)]] + 
            Sin[\((x1[t] /. aa)\)[\([1]\)] + \((x2[t] /. 
                      aa)\)[\([1]\)]], 	
	\(-	
              Cos[\((x1[t] /. aa)\)[\([1]\)]]\) - 
            Cos[\((x1[t] /. aa)\)[\([1]\)] + \((x2[t] /. 
                      aa)\)[\([1]\)]]}, {t, 0, tmax}, 
PlotPoints -> 2000, 
        Frame -> True];\)\)\)