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)')