function sol = prob2
% prob2, prob2f 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];

for tau = [1 7.5]
   sol = dde23(@prob2f,tau,history,[0, 350]);
   figure
   plot(sol.x,sol.y(1,:))
   title(['Problem 2. Baroflex Feedback Mechanism with' ...
          ' \tau = ',num2str(tau),'.'])
   xlabel('time t')
   ylabel('P_a(t)')
   axis([0 350 82 96])
end

%=============================================
function yp = prob2f(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.
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;