
function xx=Fzombie(x,alpha, beta,delta,zeta );


xx(1)=  - beta*x(1)*x(2)//-delta*x(1)
xx(2)=beta*x(1)*x(2)+zeta*x(3)-alpha*x(1)*x(2)
xx(3)=delta*x(1) + alpha*x(1)*x(2) - zeta*x(3)

endfunction

function xx=Fzombie2(x,alpha, beta,delta,zeta );

c = 0.008;
xx(1)=  - beta*x(1)*x(2)+c*x(2)
xx(2)=beta*x(1)*x(2)+zeta*x(3)-alpha*x(1)*x(2)-c*x(2)
xx(3)=delta*x(1) + alpha*x(1)*x(2) - zeta*x(3)

endfunction


function yy=eulE(y0,dt,T,alpha, beta,delta,zeta );
y=y0;
yy=y;
N=floor(T/dt);

for i=1:N,
    y=y+dt*Fzombie(y,alpha, beta,delta,zeta)
    yy=[yy,y];
end
endfunction


function yy=eulE2(y0,dt,T,alpha, beta,delta,zeta );
y=y0;
yy=y;
N=floor(T/dt);

for i=1:N,
    y=y+dt*Fzombie2(y,alpha, beta,delta,zeta)
    yy=[yy,y];
end
endfunction




