//Conditions initiales
s0=1;
r0=2;

//Fonction d'évolution des sardines
function s=sardines(t,s0,r0)
    s=2*(s0-r0)*exp(3*t)-(s0-2*r0)*exp(2*t);
endfunction

//Fonction d'évolution des requins
function r=requins(t,s0,r0)
    r=(s0-r0)*exp(3*t)-(s0-2*r0)*exp(2*t);
endfunction

//Echantillonage de t
dt=0.05
t=0:dt:1;

// Calcul par l'expression analytiques de s(t) et r(t) (méthode 1)
s=sardines(t,s0,r0);
r=requins(t,s0,r0);

// Calcul par l'exponentielle de la matrice (méthode 2)
A=[4 -2;1 1]    //Matrice du système d'équation différentielle

X=zeros(2,length(t));
for i=1:length(t)
  X(:,i)=expm(A*t(i))*[s0;r0];
end

//Tracé des courbes (t, s(t)) et (t,r(t))
xset("window", 0);
xtitle ('Evolution des sardines et des requins en fn du temps','t','s(t) et r(t)');
plot(t,s,'r')
plot(t,r,'b')

//Tracé de la courbe (r(t), s(t))
xset("window", 1);
xtitle ('Evolution relative des sardines et des requins','s(t)','r(t)');
plot(s,r)

xset("window", 2);
xtitle ('Evolution des sardines et des requins en fn du temps, calcul par exponentielle de la matrice','s(t) et r(t)');
plot(t,X(1,:),'r')
plot(t,X(2,:),'b')

xset("window", 3);
xtitle ('Evolution relative des sardines et des requins, calcul par exponentielle de la matrice','s(t)','r(t)');
plot(X(1,:),X(2,:));


