
delayp = 1.7; 
delaym = 7.1;
delays = [delayp delaym];

eps = 0.2;
ti = 0 ;
tf = 500;

T = delayp+delaym; 
c = 0.23; 
pcrit = 40; 
k = 33; 
a = 4.5;
b = 0.23;

ndelaym = floor(delaym/eps);
ndelayp = floor(delayp/eps);
nfinal = floor(tf/eps);          

% Initial conditions 
m(1)=0;
p(1)=0;    
tp(1)=0;

for i=2:2500
    
    if (i > ndelayp) 
        p(i) = p(i-1) + eps*(a*m(i-ndelayp) - b*p(i-1));
    else
        p(i) = p(i-1) + eps*(-b*p(i-1));
    end
    
    if (i > ndelaym) 
        m(i) = m(i-1) + eps*(k/(1 + (p(i-ndelaym)/pcrit)^2) - c*m(i-1));
      
    else
        m(i) = m(i-1) + eps*(k -c*m(i-1));     
        
    end
  
    tp(i) = i*eps; 
    
end 
p(1:20)
figure;

% Add code here to plot figure 

