
function [t,y]=solution(a1,a2,thetaInit,uInit,tmax,te),
    omega1=2*%pi*a1;
    omega2=2*%pi*a2;
    function deriv=systeme(t,y),
        deriv(1)=y(2),
        deriv(2)=-omega1^2*sin(y(1))-omega2^2*sin(y(1)-2*%pi*t),
    endfunction
    function J=jacobien(t,y), 
        J(1,1)=0,
        J(1,2)=1, 
        J(2,1)=-omega1^2*cos(y(1))-omega2^2*cos(y(1)-2*%pi*t), 
        J(2,2)=0,
    endfunction  
    t=[0:te:tmax];
    tolA=1d-7;
    tolR=1d-10;
    y=ode([thetaInit;uInit],0,t,tolR,tolA,systeme,jacobien);
endfunction
             
function poincare(a1,a2,initial,tmax,rect),
    s=size(initial);
    for k=1:s(1),
        [t,y]=solution(a1,a2,initial(k,1),initial(k,2),tmax,1);
        plot2d(y(1,:),y(2,:),style=0,rect=rect);
    end
endfunction
            
function spectre(a1,a2,thetaInit,uInit,tmax,te,dbmin,dbmax),
    [t,y]=solution(a1,a2,thetaInit,uInit,tmax,te);
    tfd=fft(y(1,:));
    n=length(y(1,:));
    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,dbmin,1/te/2,dbmax]); 
endfunction
            