ex03bch4.m

these codes are for solving OED with matlab

function sol = ex03bch4 
options = ddeset('RelTol',1e-5);
for tau = [7, 9]
    sol = dde23(@ddes,tau,[44; 2],[0,250],options);
    figure
    plot(sol.y(1,:),sol.y(2,:),40,12,'r*');
    title(['Plot in phase plane for \tau = ',num2str(tau),'.'])
    xlabel('Equilibrium point at (40,12).')
end
%============================================================
function dydt = ddes(t,y,Z)
xoft = y(1);
yoft = y(2);
xlag = Z(1);
dydt = [ xoft*(2*(1 - xoft/50) - yoft/(xoft + 40)) - 10
         yoft*(-3 + 6*xlag/(xlag + 40))                 ];