alpha = 0.005;
beta = 0.0095;
zeta = 0.0001;
delta = 0.0001;

N = 500;
y0(1) = N;
y0(2) = 0;
y0(3) = 0;
T = 30;
dt = 0.01;

// premier modele
yy=eulE(y0,dt,T);

// deuxieme modele
//yy=eulE2(y0,dt,T);
time = [0:dt:dt*floor(T/dt)];
plot(time',yy(1,:)');


plot2d(time',[yy(1,:)',yy(2,:)',yy(3,:)'],[1,2,3],leg="S@Z@D)");


