ex03ach4.m

these codes are for solving OED with matlab

function ex03ach4 
tspan = [0,100];
options = ddeset('RelTol',1e-5);
sol = dde23(@odes,[],[40; 16],tspan,options);
plot(sol.y(1,:),sol.y(2,:),20,15,'r*');
title('Plot in phase plane for \tau = 0.')
xlabel('Equilibrium point at (20,15).')
sol = dde23(@ddes,0.826,[40; 2],tspan,options);
figure
plot(sol.y(1,:),sol.y(2,:),20,15,'r*');
title('Plot in phase plane for \tau = 0.826')
xlabel('Equilibrium point at (20,15).')
%==========================================================
% Functions ODES and DDES differ only in the lag term XLAG.
function dydt = odes(t,y,Z)
xoft = y(1);
yoft = y(2);
xlag = xoft;
dydt = [ xoft*(2*(1 - xlag/40) - yoft/(xoft + 10)) - 10
         yoft*(xoft/(xoft + 10) - 2/3)                  ];

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