  
solution[thetaInit_,tmax_]:=Module[{sys}, 
    sys={theta'[t]==dtheta[t],dtheta'[t]==-4*Pi^2*Sin[theta[t]],
    theta[0]==thetaInit,dtheta[0]==0};
    NDSolve[sys,{theta,dtheta},{t,0,tmax},
            Method->{"ExplicitRungeKutta",DifferenceOrder->4},
            AccuracyGoal->5,MaxSteps->Infinity]
]
            

spectre[solution_,tmax_,te_]:=Module[{n,echant,tfd,sp},
    n=tmax/te;
    echant=Table[First[theta[(k-1)*te]/.solution],{k,n}];
    tfd=Fourier[echant,FourierParameters->{-1,-1}];
    sp=Table[{(k-1)/tmax,20*Log[10,Abs[tfd[[k]]]]},{k,n}];
    ListPlot[sp,PlotRange->{{0,1/te/2},{-100,0}},Joined->True]
] 
            
tmax=100;
s1=solution[0.1,tmax]
p1=ParametricPlot[{theta[t],dtheta[t]/(2*Pi)}/.s1,{t,0,10},AxesLabel->{"theta","dtheta"}]
spectre[s1,tmax,0.1]
s2=solution[Pi/2,tmax]
p2=ParametricPlot[{theta[t],dtheta[t]/(2*Pi)}/.s2,{t,0,10},AxesLabel->{"theta","dtheta"}]
spectre[s2,tmax,0.1]
s3=solution[3.1,tmax]
p3=ParametricPlot[{theta[t],dtheta[t]/(2*Pi)}/.s3,{t,0,10},AxesLabel->{"theta","dtheta"}]
spectre[s3,tmax,0.1]
Plot[0.5*dtheta[t]^2-(2*Pi)^2*Cos[theta[t]]/.s3,{t,0,tmax},PlotRange->{{0,tmax},{30,50}},
                    AxesLabel->{"t","Energie"}]
Show[p1,p2,p3,PlotRange->{{-3,3},{-3,3}}]
