r1=1; K1=3;


a=1; b=1;
m=1;

// Lokta voltera

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

//avec terme de décroissance pour les prédateurs et logistique pour les proies
// si m/b > K1  un pt d'equilibre stable (K,0), ie plus de prédateurs!
// si m/b <K1   un pt stable (m/b, r1/a(1-m/K1/b))

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

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

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,:));
plot2d(m/b,r1/a*(1-m/K1/b),-3)
title("diagramme de phase (proie,predateurs)");





// avec terme logistique pour les predateurs aussi

r2=1; K2=3;
a=0.05; b=0.1;

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


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

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


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

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

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


