function sol = prob_4_13
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% First model:
%
opts = ddeset( 'RelTol', 1e-5 );
sol1 = dde23(@odes,[],[40 16],[0,100],opts);
figure; plot( sol1.y(1,:), sol1.y(2,:), '-og' ); grid on;
sol2 = dde23(@ddes1,[0.826],[40 2],[0,100],opts);
figure; plot( sol2.y(1,:), sol2.y(2,:), '-og' ); grid on;

% Second model:
%
opts = ddeset( 'RelTol', 1e-5 );
sol1 = dde23(@ddes2,[7],[44 2],[0,250],opts);
figure; plot( sol1.y(1,:), sol1.y(2,:), '-og' ); grid on;
sol2 = dde23(@ddes2,[9],[42 2],[0,250],opts);
figure; plot( sol2.y(1,:), sol2.y(2,:), '-og' ); grid on;

% Third model:
%
opts = ddeset( 'RelTol', 1e-5 );
sol1 = dde23(@ddes3,[0.05],[2 10],[0,10],opts);
figure; plot( sol1.y(1,:), sol1.y(2,:), '-og' ); grid on;

return;

function dydt = odes(t,y,Z)

dydt = [ y(1)*( 2*(1-y(1)/40) - y(2)/(y(1)+10) ) - 10;
         y(2)*( y(1)/(y(1)+10) - 2/3 );        
         ];

function dydt = ddes1(t,y,Z)

dydt = [ y(1)*( 2*(1-Z(1,1)/40) - y(2)/(y(1)+10) ) - 10;
         y(2)*( y(1)/(y(1)+10) - 2/3 );        
         ];

function dydt = ddes2(t,y,Z)

dydt = [ y(1)*( 2*(1-y(1)/50) - y(2)/(y(1)+40) ) - 10;
         y(2)*( -3 + 6*Z(1,1)/(Z(1,1)+40) );
         ];

function dydt = ddes3(t,y,Z)

dydt = [ y(1)*( 20 - y(1) - y(2) ) - 7;
         -15*y(2) + 3*Z(1,1)*Z(2,1);
         ];