ex09ch4.m

these codes are for solving OED with matlab

function sol = ex09ch4
global h1 h2 h3 h4 h5 h6 h7 h8
h1 = 2; h2 = 0.8; h3 = 1e4; h4 = 0.17;
h5 = 0.5; h6 = 300; h7 = 0.12; h8 = 8;
state = +1;
opts = ddeset('Jumps',-1e-6,'Events',@events,...
              'RelTol',1e-5,'AbsTol',1e-8);
sol = dde23(@ddes,0.5,@history,[0, 60],opts,state);
while sol.x(end) < 60
   state = -state;
   sol = dde23(@ddes,0.5,sol,[sol.x(end), 60],opts,state);
end
yplot = sol.y;
yplot(1,:) = 1e4*yplot(1,:);
yplot(2,:) = yplot(2,:)/2;
yplot(4,:) = 10*yplot(4,:);
plot(sol.x,yplot)
legend('1e4*y_1','y_2/2','y_3','10*y_4',0)
axis([0 60 -1 15.5])
%=========================================================
function dydt = ddes(t,y,Z,state)
global h1 h2 h3 h4 h5 h6 h7 h8
if state == +1
   xi = 1;
else
   xi = (1 - y(4))*10/9;
end
ylag = Z;
dydt = zeros(4,1);
dydt(1) = (h1 - h2*y(3))*y(1);
dydt(2) = xi*h3*ylag(3)*ylag(1) - h5*(y(2) - 1);
dydt(3) = h4*(y(2) - y(3)) - h8*y(3)*y(1);
dydt(4) = h6*y(1) - h7*y(4);

function v = history(t,state)
v = [1; 1; 1; 0];
v(1) = max(0,1e-6 + t);

function [value,isterminal,direction] = events(t,y,Z,state)
value = y(4) - 0.1;
isterminal = 1;
direction = 0;