function sol = prob_4_5(G)
% 
% Written by:
% -- 
% John L. Weatherwax                2005-04-27
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

global c D n
c = 0.5; D = 5; G = 1; n = 100; 
% x = y(1), y = y(2), lambda = y(3), 
% I_1 = y(4), I_2 = y(5), I_3 = y(6).
y0 = [0.8*n; 0.2*n; 0; 0; 0; 0];
sol = dde23(@odes,[],y0,[0, D],[],G);
sol = dde23(@ddes,D,sol,[D, 4*D],[],G);
plot(sol.x,[sol.y(1:2,:); 100*sol.y(3,:)]);
legend('x(t)','y(t)','100 \lambda(t)',0)
disp( 'Equilibrium solution #1:' );
disp( [ n 0 0 ] );
disp( 'Equilibrium solution #2:' );
disp( [ G*n/(c*D) n*(1-G/(c*D)) c*D-G ] );
disp( 'The long time solution limit''s is:' );
sol.y(1:3,end)
%============================================================
function dydt = odes(t,y,Z,G)
global c D n
dydt = zeros(6,1);
dydt(1) = - y(1)*y(3) + G*y(2);
dydt(2) = - dydt(1);
dydt(4) = exp(G*t)*y(2);
dydt(5) = t*exp(G*t)*y(1)*y(3);
dydt(6) = exp(G*t)*y(1)*y(3);
dydt(3) = (c/n)*exp(-G*t)*((dydt(4)+dydt(5))-G*(y(4)+y(5)));

function dydt = ddes(t,y,Z,G)
global c D n
dydt = zeros(6,1);
dydt(1) = - y(1)*y(3) + G*y(2);
dydt(2) = - dydt(1);
dydt(4) = exp(G*t)*y(2) - exp(G*(t - D))*Z(2);
dydt(5) = D*exp(G*t)*y(1)*y(3) - y(6);
dydt(6) = exp(G*t)*y(1)*y(3) - exp(G*(t - D))*Z(1)*Z(3);
dydt(3) = (c/n)*exp(-G*t)*((dydt(4)+dydt(5))-G*(y(4)+y(5)));