ex10bch4.m

these codes are for solving OED with matlab

function sol = ex10bch4
global r Vstr
R = 1.05; r = 0.068; Vstr = 67.9;

P0      = 93;
Paval   = P0;
Pvval   = (1 / (1 + R/r)) * P0;
Hval    = (1 / (R * Vstr)) * (1 / (1 + r/R)) * P0;
history = [Paval; Pvval; Hval];
tau = 4;
opts = ddeset('Jumps',600,'RelTol',1e-4,'AbsTol',1e-7);
sol = dde23(@ddes,tau,history,[0, 1000],opts);
plot(sol.x,sol.y(3,:))
title(['Baroflex Feedback Mechanism with' ...
       ' \tau = ',num2str(tau)])
xlabel('time t')
ylabel('H(t)')

%====================================================
function v = ddes(t,y,Z)
global r Vstr
ca = 1.55; cv = 519; gammaH = 0;
alpha0 = 93; alphas = 93; alphap = 93; alphaH = 0.84;
beta0  = 7; betas  = 7; betap  = 7; betaH  = 1.17;

if t <= 600
   R = 1.05;
else
   R = 0.21 * exp(600-t) + 0.84;
end    
ylag = Z(:,1);
Patau = ylag(1);
Paoft = y(1);
Pvoft = y(2);
Hoft  = y(3);
dPadt = - (1 / (ca * R)) * Paoft + (1/(ca * R)) * Pvoft ...
        + (1/ca) * Vstr * Hoft;
dPvdt = (1 / (cv * R)) * Paoft                          ...
        - ( 1 / (cv * R) + 1 / (cv * r) ) * Pvoft;
Ts = 1 / ( 1 + (Patau / alphas)^betas );
Tp = 1 / ( 1 + (alphap / Paoft)^betap );
dHdt = (alphaH * Ts) / (1 + gammaH * Tp) - betaH * Tp;
v = [ dPadt; dPvdt; dHdt ];