ex32ch2.m

these codes are for solving OED with matlab

function ex32ch2
% 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;
options = odeset('Mass',@mass,'MassSingular','no','Events',@events);
tD = linspace(0,1,11);
zU = zeros(11,1);
for i = 1:11
    G0 = 270.9*(0.8 + 0.2*exp(-tD(i)));
    [z,y,ze,ye,ie] = ode45(@odes,[0 5],[795.5; 255.0],options);
    zU(i) = ze(end);
end 
plot(tD,zU); 
title('Upper boundary of the liquid region.')
%===================================================================
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;