ex13bch4.m

these codes are for solving OED with matlab

function sol = ex13bch4 
% Eq. (3.11) / Figs. 3.2-3.3
% Solution oscillates around, and converges to the equilibrium point (40,12)
% for tau=7 and approaches a limit cycle for tau=9
for tau = [7 9]
options = ddeset('RelTol', 1.e-5, 'AbsTol', 1.e-5);
   sol = dde23(@ddes, tau, [44; 2], [0, 250], options);
   figure
   plot(sol.x(:), sol.y(1,:), sol.x(:), sol.y(2,:));
   legend('x(t)', 'y(t)')
   figure
   plot(sol.y(1,:), sol.y(2,:))
   title(['Phase plot for', ' \tau = ', num2str(tau), '.'])
end
%============================================================

function dydt = ddes(t,y,Z)
dydt = zeros(2,1);
xoft = y(1);
yoft = y(2);
xlag = Z(1,1);
dydt(1) = xoft*(2*(1-xoft/50)-yoft/(xoft+40))-10;
dydt(2) = yoft*(-3+(6*xlag)/(xlag+40));