function v = 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;
v = zeros(3,1);
v(1) = dPadt;
v(2) = dPvdt;
v(3) = dHdt;