r1=1; K1=3;

D=1;
a=1; b=2;
m=1;

// Holling

//avec terme de décroissance pour les prédateurs et logistique pour les proies
// si m*D/(b-m) > K1  un pt d'equilibre stable (K,0), ie plus de prédateurs!
// si m*D/(b-m) < K1   compliqué! 
// 3 cas: 1) cycle limite avec oscillations: ex: r1=1, K1=3, D=1,a=1,b=2,m=1
//  2) un etat stable attractif:  ex r1=1, K1=1.2, D=1,a=1,b=2,m=1
//  3) extenction predateurs (K,0) j'y arrive pas!

deff('[yprime]=f(t,y)',['yprime(1)=r1*y(1)*(1-y(1)/K1) - a*y(1)*y(2)/(y(1)+D)' 'yprime(2)= - m*y(2) + b*y(1)*y(2)/(y(1)+D)' ]);

y0=[1 1.]';
t0=0;
t=0:0.01:1000;

y=ode("rk",y0,t0,t,f);


xset("window",0);
clf();
plot2d(t,y(1,:));
plot2d(t,y(2,:),-1);
title("variation temporelles des proies (-) et predateurs (+) ");

xset("window",1);
clf();
plot2d(y(1,:),y(2,:));

title("diagramme de phase (proie,predateurs)");




