ex08ch4.m

these codes are for solving OED with matlab

function sol = ex08ch4
r = 0.5;
c = 1/sqrt(2);
options = ddeset('Jumps',[(1-c), 1, (2-c)],...
                 'RelTol',1e-5,'AbsTol',1e-8);
sol = dde23(@ddes,1,10,[0, 10],options,r,c);
y10 = deval(sol,10);
fprintf('DDE23 computed     y(10) =%15.11f.\n',y10);
fprintf('Reference solution y(10) =%15.11f.\n',0.06302089869);
plot(sol.x,sol.y)
title(['Hoppensteadt-Waltman model with r = ',...
       num2str(r),'.'])
xlabel('time t')
ylabel('y(t)')
Ioft = -(1/r)*(sol.yp ./ sol.y);
figure
plot(sol.x,Ioft)
title(['Hoppensteadt-Waltman model with r = ',...
       num2str(r),'.'])
xlabel('time t')
ylabel('I(t)')
%============================================================
function v = ddes(t,y,Z,r,c)
mu = r/10;
if     t <= 1 - c
   v = -r*y*0.4*(1 - t);
elseif t <= 1
   v = -r*y*(0.4*(1 - t) + 10 - exp(mu)*y);
elseif t <= 2 - c
   v = -r*y*(10 - exp(mu)*y);
else
   v = -r*exp(mu)*y*(Z - y);
end