ex13cch4.m

these codes are for solving OED with matlab

function sol = ex13cch4 
% Eq. (4.15) / Fig. 4.2
options = ddeset('RelTol', 1.e-5, 'AbsTol', 1.e-5);
sol = dde23(@ddes, 0.05, [2; 10], [0, 10], 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 = 0.05', '.'])
%============================================================

function dydt = ddes(t,y,Z)
dydt = zeros(2,1);
xoft = y(1);
yoft = y(2);
xlag = Z(1,1);
ylag = Z(2,1);
dydt(1) = xoft*(20-xoft-yoft)-7;
dydt(2) = -15*yoft+3*xlag*ylag;