k1 = 5
k2 = 10
m1 = 0.1
m2 = 0.1
L1 = 6
L2 = 6
A = [(k1+k2)/m1 -k2/m1; -k2/m2 k2/m2];

function [lambda1, p1] = puissancesIterees(A, x0, tol)
    xk = x0;
    while norm(A*xk - (xk'*A*xk)*xk)>tol
       xk = (A*xk) / norm(A*xk);
    end; // while
    lambda1 = xk'*A*xk;
    p1 = xk;
endfunction

x0 = [1;0];
tol = 1e-6;
[lambda1, p1] = puissancesIterees(A, x0, tol)
disp(lambda1);
disp(p1);

function [lambda2, p2] = puissancesInverses(A, x0, tol)
    xk = x0;
    zk = A\xk;
    // Test ||Ax-lambda*x||>tol ?
    while norm(A*xk - (xk'*A*xk)*xk)>tol
       zk = A\xk; 
       xk = (zk) / norm(zk);
    end; // while
    // xk est à la fois vecteur propre de A et A^{-1}
    // A xk = lambda2 xk
    // donc lambda2 = xk'*A*xk
    lambda2 = xk'*A*xk;
    p2 = xk;
endfunction

x0 = [1;0];
tol = 1e-6;
[lambda2, p2] = puissancesInverses(A, x0, tol)
disp(lambda2);
disp(p2);

// Verif
[V, D] = spec(A)
disp(V), disp(D)