function sol = prob_4_11
% 
% 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.
% 
%-----

global a b r;
a = 0.8; b = 0.7; r = 0.08;

% Compute y_1_0 and y_2_0:
global y_1_0 y_2_0; 
% Steady state solution for y_2_0:
C = [ -((b^3)/3) (b^2)*a (b-1-b*(a^2)) (-a + (a^3)/3) ];
sol = roots( C ); 
y_2_0 = sol(3);
y_1_0 = b*y_2_0 - a;
y_1_0 = -1.22764016;

tau = 20; m = +10;
sol = dde23(@ddes,[tau], @history , [0, 25], [], m );
figure; plot( sol.x, sol.y(1,:), '-o' ); 

tau = 20; m = -10;
sol = dde23(@ddes,[tau], @history , [0, 60], [], m );
figure; plot( sol.x, sol.y(1,:), '-o' ); 




function v = history(t,m)

global y_1_0 y_2_0;  
  
v = [ 0.4*y_1_0; 1.8*y_2_0 ];


function dydt = ddes(t,y,Z,m)

global a b r y_1_0 y_2_0;
dydt = [ y(1) - (y(1)^3)/3 - y(2) + m*( Z(1,1) - y_1_0 );
         r*( y(1) + a - b*y(2) );
         ];