s0=10000; // Conditions initiales
r0=200;

dt=0.005;
t=0:dt:4;
n=length(t);
a=2.1;b=10.2;c=0.0013;d=0.0024;

// 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(1,k+1)=Y(1,k)+dt*(a*Y(1,k) - c*Y(1,k)*Y(2,k));
    Y(2,k+1)=Y(2,k)+dt*(-b*Y(2,k) + d*Y(1,k)*Y(2,k));
end

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

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

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