ex10ach4.m

these codes are for solving OED with matlab

function sol = ex10ach4
global R 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];
opts = ddeset('RelTol',1e-5,'AbsTol',1e-8);
for tau = [1 7.5]
   sol = dde23(@ddes,tau,history,[0, 350],opts);
   figure
   plot(sol.x,sol.y(1,:))
   title(['Baroflex Feedback Mechanism with' ...
          ' \tau = ',num2str(tau)])
   xlabel('time t')
   ylabel('P_a(t)')
   axis([0 350 70 105])
end
%=======================================================
function v = ddes(t,y,Z)
global R 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;
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 ];