function sol = prob_4_10
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

global c_a c_v R r alpha_0 alpha_s alpha_p alpha_H beta_0 beta_s beta_p beta_H gamma_H V_str P_0;
c_a = 1.55; 
c_v = 519;
R = 1.05;
r = 0.068;
alpha_0 = 93;
alpha_s = 93;
alpha_p = 93;
alpha_H = 0.84;
beta_0 = 7;
beta_s = 7;
beta_p = 7;
beta_H = 1.17;
gamma_H = 0;
V_str = 67.9;
P_0 = 93;

options = ddeset( 'Jumps', [600] );

tau = 1; 
sol = dde23(@ddes,[tau], @history , [0, 350], options );
figure; plot( sol.x, sol.y(1,:), '-o' ); 

tau = 7.5; 
sol = dde23(@ddes,[tau], @history , [0, 350], options );
figure; plot( sol.x, sol.y(1,:), '-o' ); 

tau = 4; 
sol = dde23(@ddes,[tau], @history , [0, 1000], options );
figure; plot( sol.x, sol.y(1,:), '-o' ); 


function v = history(t)

global c_a c_v R r alpha_0 alpha_s alpha_p alpha_H beta_0 beta_s beta_p beta_H gamma_H V_str P_0;

v = [ P_0;
     ( r/(r+R) )*P_0;
     ( 1/(r+R) )*(P_0/V_str);
    ];


function dydt = ddes(t,y,Z)

global c_a c_v R r alpha_0 alpha_s alpha_p alpha_H beta_0 beta_s beta_p beta_H gamma_H V_str P_0;

Ts = ( 1 + (Z(1,1)/alpha_s)^beta_s )^(-1);
Tp = ( 1 + (alpha_p/y(1))^beta_p )^(-1);
f_Ts_Tp = ( alpha_H*Ts )/( 1 + gamma_H*Tp ) - beta_H*Tp;

if( t <= 600 ) 
  R = 1.05;
else
  R = 0.21 * exp( 600-t ) + 0.84;
end

dydt = [ -1/(c_a*R)*y(1) + (1/(c_a*R))*y(2) + (V_str/c_a)*y(3); 
         (1/(c_v*R))*y(1) - ( 1/(c_v*R) + 1/(c_v*r) )*y(2);
         f_Ts_Tp;
         ];