ex13ach4.m

these codes are for solving OED with matlab

function sol = ex13ach4 
% Eq. (2.18) / Figs. 2.1-2.2
% Solution approaches equilibrium point (20,15)
options = ddeset('RelTol' ,1.e-5, 'AbsTol', 1.e-5);
sol = dde23(@odes, [], [40; 16], [0, 100], 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', '.'])

% Eq. (2.18) / Fig. 2.4
% Solution oscillates periodically around the equilibrium point (20,15)
options = ddeset('RelTol', 1.e-5, 'AbsTol', 1.e-5);
sol = dde23(@ddes, 0.826, [40 2], [0, 100], 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.826', '.'])
%============================================================

function dydt = odes(t,y,Z)
global r a b H J K
dydt = zeros(2,1);
xoft = y(1);
yoft = y(2);
xlag = xoft;
dydt(1) = xoft*(2*(1-xlag/40)-yoft/(xoft+10))-10;
dydt(2) = yoft*(xoft/(xoft+10)-2/3);

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-xlag/40)-yoft/(xoft+10))-10;
dydt(2) = yoft*(xoft/(xoft+10)-2/3);