z0=1;
y0=0.5;

dt=0.001;
epsl=0.001;

N=round(1/dt);

z=zeros(N+1,1);
y=zeros(N+1,1);
t=zeros(N+1,1);

z(1)=z0;
y(1)=y0;
t(1)=0;

for n=1:N
    t(n+1)=n*dt;
    y(n+1)=(y(n)+dt/epsl)/(1+dt/epsl);
    z(n+1)=z(n)/(1-dt*y(n+1));
end

plot(t,y)
hold on
plot(t,z)
% axis([-0.1 1.1 -0.1 1.1])
legend('y','z')