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

A=[4 -2;1 1]

dt=0.005;
t=0:dt:10;
n=length(t);
a=3;b=1;c=2;d=1;

// Solution approximative méthode d'Euler (1)
Y=zeros(2,n);    // préparation de la matrice (1ère L : sardines, 2ème L : requins)
Y(1,1)=s0;    // condition initiale sardines
Y(2,1)=r0;    // condition initiale requins

// remplissage de la matrice (méthode d'Euler)
for k=1:n-1
    Y(:,k+1)=Y(:,k)+dt*A*Y(:,k);
end

// Solution analytique par calcul de l'exponentielle de la matrice (2)
X=zeros(2,length(t));
for i=1:length(t)
  X(:,i)=expm(A*t(i))*[s0;r0];
end

//Les sardines
xset("window", 0);
xtitle ('Evolution des sardines analytique et numérique','t','s(t) et s_num(t)');
plot(t,Y(1,:),'r')
plot(t,X(1,:),'b')

//Les requins
xset("window", 1);
xtitle ('Evolution des requins analytique et numérique','t','r(t) et r_num(t)');
plot(t,Y(2,:),'r')
plot(t,X(2,:),'b')
