p=1.5;
a=0.1;
b=8;
dj=1.2;
tau=1;
d=0.6;
m=0.1;

history = [0.5; 1];
tspan = [0, 100]; 
opts = ddeset('RelTol',1e-5,'AbsTol',1e-8);

% Solve equation
sol = dde23('stagepp', tau, history, tspan, opts,p,a,b,dj,tau,d,m);

plot(sol.x, sol.y);
legend('prey', 'predator', 1);
title('Figure 1. A predator-prey model with stage structure');
text(10,1.93,'x\prime=x(1-x)-ypx/(1+ax)','FontSize',10)
text(10,1.83,'y\prime=be^{-d_j\tau}py(t-\tau)x(t-\tau)/(1+ax(t-\tau))-dy-my^2','FontSize',10)
text(10,1.7,'p=1.5, a=0.1, b=8, d_j=1.2, \tau=1, d=0.6, m=0.1','FontSize',10)
text(10,1.6,'(x(\theta), y(\theta))=(0.5, 1), \theta \in [-\tau, 0]','FontSize',10)
ylim([0, 2]);
xlabel('time t');
ylabel('x, y');