function prob_2_32
% 
% 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.
% 
%-----
  
% Define the physical constants:
global kappa beta a Cp K ga sinth Phi Ph Af G0
kappa = 0.171446272015689e-8;
beta  = 0.213024626664637e-2;
a     = 0.108595374561510e+4;
Cp    = 0.496941623289027e+4;
K     = 10;
ga    = 9.80665;
sinth = 1;
Phi   = 1.1e+5;
Ph    = 797.318;
Af    = 3.82760;
%G0    = 270.9;
tD = linspace( 0, 1, 11 );
G0Array = 270.9*( 0.8 + 0.2*exp(-tD) );

options = odeset('Mass',@mass,'MassSingular','no','Events',@events);

zU = [];
for g0i=G0Array,
  G0 = g0i;
  [z,y,ze,ye,ie] = ode45(@odes,[0 5],[795.5; 255.0],options);  
  zU = [zU, ze(end)]; 
end
plot(tD,zU,'-o');
%===================================================================
function dydz = odes(z,y)
global kappa beta a Cp K ga sinth Phi Ph Af G0
rho = y(1);
T = y(2);
dydz = [ (-K*G0*abs(G0/rho) - rho*ga*sinth)
         (a^2 *Phi*Ph*kappa)/(Cp*Af)       ];
 
function A = mass(z,y)
global kappa beta a Cp K ga sinth Phi Ph Af G0
rho = y(1);
T = y(2);
A = zeros(2);
A(1,1) = 1/(rho*kappa) - (G0/rho)^2;
A(1,2) = beta/kappa;
A(2,1) = -(a^2 *beta*(T + 273.15)*G0)/(Cp*rho^2);
A(2,2) = G0/rho;

function [value,isterminal,direction] = events(z,y)
isterminal = 1;
direction  = 0;
rho   = y(1);
T     = y(2);
rhosat  = -3.3*(T - 290.0) + 738.0;
value = rho - rhosat;