b = 2;
c = 1;
e = 1 - c/b;   % Equilibrium point.

sol = dde23('prob4f',7,0.8,[0, 60],[],b,c);

plot(sol.x,sol.y,[0, 60],[e, e],'r')
legend('y(t)','1 - c/b')
title(['Problem 4. Cooke Epidemic Model with' ...
       ' b = ',num2str(b),' , c = ',num2str(c),'.'])
xlabel('time t')