// Définition d'une fonction construisant Qi+1
function[Qip1]=DefQip1(i,u,y,N)
    Qip1=zeros(2*(i+1),2*(i+1));// déclare Qip1 comme une matrice 
                                // 2(i+1)x2(i+1) remplie de zéros
    for k=i+1:N-(i+1),
        A=zeros(2*(i+1),1); // déclare le vecteur col. dans Qi+1
        B=zeros(1,2*(i+1)); // déclare le vecteur lig. dans Qi+1
        // construction de A et B
        for j=1:i+1, 
            A(2*j-1,1)= u(1,k-j+1);//lignes d'indices impairs dans A
            A(2*j,1)= u(1,k+j);//lignes d'indices pairs dans A
            B(1,2*j-1)= y(1,k+j);//col. d'indices impairs de B
            B(1,2*j)= u(1,k+j);//col. d'indices pairs de B
         end 
         Qip1=Qip1+A*B;      
    end
    Qip1=Qip1/N;
endfunction

load("Qu1_1.sod");
N=size(u,2); // affecte à N le nbre de col. de u
// On construit Qi pour i=1 en appelant DefQip1 pour i=0
Qi=DefQip1(0,u,y,N);
for i=1:15; // On pense que l'ordre du système est inf. à 15
    // Construction de Qi+1
    Qip1=DefQip1(i,u,y,N);
    // Calcul et affichage dy RDI
    disp(i)
    disp(abs(det(Qi)/det(Qip1)));
    Qi=Qip1;
end