(* Double Pendulum Numerical integration *)
(* θ is the coordinate on the top pendulum, φ is the coordinate on the lower pendulum. Each pendulum has the same length *)

x[θ_,φ_]:= 2-(Cos[θ-φ])^2

eqn5 := -(((pθ[t])^2)/2 + (pφ[t])^2 - Cos[θ[t]-φ[t]]*pθ[t]*pφ[t])*(2*Cos[θ[t]-φ[t]]*Sin[θ[t]-φ[t]])/((x[θ[t],φ[t]])^2) -Sin[θ[t]-φ[t]]*pθ[t]*pφ[t]/x[θ[t],φ[t]] -2Sin[θ[t]]

eqn6 := (pθ[t]-Cos[θ[t]-φ[t]]*pφ[t])/x[θ[t],φ[t]]

eqn7 := -((pθ[t]^2)/2 + pφ[t]^2 - Cos[θ[t]-φ[t]]*pθ[t]*pφ[t])*(2*Cos[θ[t]-φ[t]]*Sin[θ[t]-φ[t]])/(x[θ[t],φ[t]]^2) +Sin[θ[t]-φ[t]]*pθ[t]*pφ[t]/x[θ[t],φ[t]] -Sin[φ[t]]

eqn8 := (2*pφ[t]-Cos[θ[t]-φ[t]]*pθ[t])/x[θ[t],φ[t]]

(* Generating a solution for the initial conditions θ[0]==Pi/2 *)
sol= NDSolve[{pθ'[t]==eqn5,pφ'[t]==eqn7,θ'[t]==eqn6,φ'[t]==eqn8, pθ[0]== pφ[0]==0,φ[0]==0,θ[0]==Pi/2},{θ,pθ,φ,pφ}, {t,100},MaxSteps->Infinity]

(* Viewing the motion *)
Plot[Evaluate[φ[t] /. sol], {t,0,50}]

Plot[Evaluate[θ[t] /. sol],{t,0,50}]

(* The phase space plot for the same motion *)
ParametricPlot[Evaluate[{φ[t],pφ[t]} /. sol], {t,0,50}, PlotRange -> All]

(* To produce a pair of solutions we need expressions for the two solutions
{θ1,pθ1,φ1,pφ1,θ2,pθ2,φ2,pφ2} *) 
eqn51 := -(((pθ1[t])^2)/2 + (pφ1[t])^2 - Cos[θ1[t]-φ1[t]]*pθ1[t]*pφ1[t])*(2*Cos[θ1[t]-φ1[t]]*Sin[θ1[t]-φ1[t]])/((x[θ1[t],φ1[t]])^2) -Sin[θ1[t]-φ1[t]]*pθ1[t]*pφ1[t]/x[θ1[t],φ1[t]] -2Sin[θ1[t]]

eqn61 := (pθ1[t]-Cos[θ1[t]-φ1[t]]*pφ1[t])/x[θ1[t],φ1[t]]

eqn71 := -((pθ1[t]^2)/2 + pφ1[t]^2 - Cos[θ1[t]-φ1[t]]*pθ1[t]*pφ1[t])*(2*Cos[θ1[t]-φ1[t]]*Sin[θ1[t]-φ1[t]])/(x[θ1[t],φ1[t]]^2) +Sin[θ1[t]-φ1[t]]*pθ1[t]*pφ1[t]/x[θ1[t],φ1[t]] -Sin[φ1[t]]

eqn81 := (2*pφ1[t]-Cos[θ1[t]-φ1[t]]*pθ1[t])/x[θ1[t],φ1[t]]

eqn52 := -(((pθ2[t])^2)/2 + (pφ2[t])^2 - Cos[θ2[t]-φ2[t]]*pθ2[t]*pφ2[t])*(2*Cos[θ2[t]-φ2[t]]*Sin[θ2[t]-φ2[t]])/((x[θ2[t],φ2[t]])^2) -Sin[θ2[t]-φ2[t]]*pθ2[t]*pφ2[t]/x[θ2[t],φ2[t]] -2Sin[θ2[t]]

eqn62 := (pθ2[t]-Cos[θ2[t]-φ2[t]]*pφ2[t])/x[θ2[t],φ2[t]]

eqn72 := -((pθ2[t]^2)/2 + pφ2[t]^2 - Cos[θ2[t]-φ2[t]]*pθ2[t]*pφ2[t])*(2*Cos[θ2[t]-φ2[t]]*Sin[θ2[t]-φ2[t]])/(x[θ2[t],φ2[t]]^2) +Sin[θ2[t]-φ2[t]]*pθ2[t]*pφ2[t]/x[θ2[t],φ2[t]] -Sin[φ2[t]]

eqn82 := (2*pφ2[t]-Cos[θ2[t]-φ2[t]]*pθ2[t])/x[θ2[t],φ2[t]]

(* Generating a pair of solutions for the initial conditions φ1[0]=π and φ1[0]=π+0.001 and θ1[0]=π/2 or 90 degrees and θ2[0]=π/2 *)
sol121= NDSolve[{pθ1'[t]==eqn51,pφ1'[t]==eqn71,θ1'[t]==eqn61,φ1'[t]==eqn81, pθ1[0]== pφ1[0]==0,φ1[0]==Pi,θ1[0]== Pi/2,pθ2'[t]==eqn52,pφ2'[t]==eqn72,θ2'[t]==eqn62,φ2'[t]==eqn82, pθ2[0]== pφ2[0]==0,φ2[0]==Pi+0.001,θ2[0]== Pi/2},{θ1,pθ1,φ1,pφ1,θ2,pθ2,φ2,pφ2}, {t,1000},MaxSteps->Infinity]

Plot[Evaluate[φ1[t] /. sol121], {t,50,150}]

ParametricPlot[Evaluate[{φ1[t],pφ1[t]} /. sol121], {t,50,200}, PlotRange -> All]

Plot[Evaluate[Log[Abs[φ1[t]-φ2[t]]] /. sol121], {t,0,100}, PlotPoints->10000,PlotRange->{-18,0}]