function Vprime=F(t,V,lambda)
    w = sqrt(lambda);
    u = V(1:2,:);
    uprime = V(3:4);
    Vprime = [uprime
              -A*u+[0; 1/m2]*sin(w*t)];
endfunction

cd(get_absolute_file_path());
exec("anime_os.sci",-1)

k1=5;
k2=10;
m1=0.1;
m2=0.1;

A=[(k1+k2)/m1 -k2/m1
   -k2/m2 k2/m2];

ITMAX = 100;
TOL = 1e-8;
// puissance
x=rand(2,1);
for k = 1:ITMAX
    y = A*x;
    mu = x'*y;
    r = y-mu*x;
    if norm(r) <= TOL
        break
    end
    x = y/norm(y);
end
lambda1 = mu;

// puissance inverse
x=rand(2,1);
for k = 1:ITMAX
    y = A\x;
    mu = x'*y;
    r = y-mu*x;
    if norm(r) <= TOL
        break
    end
    x = y/norm(y);
end

lambda2= 1/mu;

[t,V] = cvode(#(t,y)->(F(t,y,lambda1)),linspace(0,10,200),zeros(4,1));
u = V(1:2,:);
clf
plot(t,u)

anime_os(t,V,m2,m2,4,4)

[t,V] = cvode(#(t,y)->(F(t,y,lambda2)),linspace(0,10,200),zeros(4,1));
u = V(1:2,:);
clf
plot(t,u)

anime_os(t,V,m2,m2,4,4)