function sol = prob2b
% prob2b, prob2bf of tutorial.

global ca cv R r Vstr gammaH       ...
       alpha0 alphas alphap alphaH ...
       beta0  betas  betap  betaH

ca     = 1.55;
cv     = 519;
R      = 1.05;
r      = 0.068;
Vstr   = 67.9;
alpha0 = 93;
alphas = 93;
alphap = 93;
alphaH = 0.84;
beta0  = 7;
betas  = 7;
betap  = 7;
betaH  = 1.17;
gammaH = 0;

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);
sol = dde23(@prob2bf,tau,history,[0, 1000],opts);
figure
plot(sol.x,sol.y(1,:))
title(['Problem 2b. Baroflex Feedback Mechanism.' ...
       ' \tau = ',num2str(tau),'.'])
xlabel('time t')
ylabel('P_a(t)')
figure
plot(sol.x,sol.y(2,:))
legend('P_v(t)')
title(['Problem 2b. Baroflex Feedback Mechanism.' ...
       ' \tau = ',num2str(tau),'.'])
xlabel('time t')
ylabel('P_v(t)')
figure
plot(sol.x,sol.y(3,:))
title(['Problem 2b. Baroflex Feedback Mechanism.' ...
       ' \tau = ',num2str(tau),'.'])
xlabel('time t')
ylabel('H(t)')

%==================================================
function yp = prob2bf(t,y,Z)
global ca cv R r Vstr gammaH       ...
       alpha0 alphas alphap alphaH ...
       beta0  betas  betap  betaH
       
% Local variables are used to express the equations in terms
% of the physical quantities of the model. 
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;
yp = zeros(3,1);
yp(1) = dPadt;
yp(2) = dPvdt;
yp(3) = dHdt;