
function [t,y,e]=solution(thetaInit,tmax,te),
    function [deriv]=systeme(t,y),
        deriv(1)=y(2),
        deriv(2)=-4*%pi^2*sin(y(1)),
    endfunction
    t=[0:te:tmax];
    tolA=1d-7;
    tolR=1d-10;
    y=ode("adams",[thetaInit;0],0,t,tolR,tolA,systeme);
    n=length(y(1,:));
    e=zeros(1,n);
    for k=1:n,
        e(k)=0.5*y(2,k)^2-4*%pi^2*cos(y(1,k));
    end;
endfunction
            

function spectre(y,tmax,te)
    tfd=fft(y);
    n=length(y);
    sp=zeros(1,n);
    f=zeros(1,n);
    for k=1:n,
        sp(k)=20*log10(abs(tfd(k)));
        f(k)=(k-1)/tmax;
    end;
    plot2d(f,sp,rect=[0,-50,1/te/2,100]);
endfunction 
            
tmax=100;
[t,y1,e1]=solution(3.1,tmax,0.01);
plotK=scf();
plot2d(t,y1(1,:),rect=[0,-%pi,5,%pi]);
xtitle("","t","theta");
plotL=scf();
plot2d(t,e1,rect=[0,30,tmax,50]);
xtitle("","t","energie");
plotI=scf();
plot2d(y1(1,:),y1(2,:))
xtitle("","theta","dtheta");
[t,y1,e1]=solution(3.1,tmax,0.1);
plotJ=scf();
spectre(y1(1,:),tmax,0.1);
xtitle("","f","|Cn|");
