f = inline('max(x,0)');

dt=0.01;
T=10;
t=0:dt:T;
L=length(t);

X=zeros(1,L);
Y=zeros(1,L);
Z=zeros(1,L);
Zcomp=zeros(1,L);

thetay=0.5;
thetaz=0.3;

for i=1:L-1
    if t(i)>2
        X(i)=1;
    end
    Y(i+1)=Y(i)+dt*(f(X(i)-thetay)-Y(i));
    Z(i+1)=Z(i)+dt*(f(X(i)-Y(i)-thetaz)-Z(i));
    Zcomp(i+1)=Zcomp(i)+dt*(f(0.5*X(i)-thetaz)-Zcomp(i));
end

figure
subplot(3,1,1)
plot(t,X,'LineWidth',1.5)
ylim([0,1.1])
ylabel('X(t)')
subplot(3,1,2)
plot(t,Y,'LineWidth',1.5)
ylim([0,0.6])
ylabel('Y(t)')
subplot(3,1,3)
hold on
plot(t,Z,'LineWidth',1.5)
plot(t,Zcomp,'r','LineWidth',1.5)
plot(t,0.2*ones(1,L),'k-.')
plot(0:0.1:2.7, 0.1,'k-.')
ylabel('Z(t)')