function Tk = outlet(fu,ox)
% First load a set of thermodynamic parameters from E&R Appendix A, Table
% 2.4 and 4.1.  Assume air is 79% N2, 21% O2.
%                 A.1      A.2      A.3         A.4    dHcomb
Thermo(1,:) = [    .56   .27559 -1.0355e-5          0 -2219e3 ]; %Propane  
Thermo(2,:) = [ 172.02 -1.08574  4.4887e-3 -6.5539e-6 -2878e3 ]; %Butane
Thermo(3,:) = [   2.02  0.44729 -1.7174e-4          0 -3509e3 ]; %Pentane
Thermo(4,:) = [  32.83 -0.03633  1.1532e-4 -1.2194e-7       0 ]; %Oxygen
Thermo(5,:) = [  31.21 -0.01680  4.2238e-5 -3.2528e-8       0 ]; %Air
Thermo(6,:) = [  18.86  0.07937 -6.7834e-5  2.4426e-8       0 ]; %Carbon Dioxide
Thermo(7,:) = [  33.80 -0.00795 -2.8228e-5 -1.3115e-8       0 ]; %Water Vapor
% The two inputs into the program are the row numbers in the matrix that
% correspond to the fuel (propane, butane, or pentane) and the oxidant
% (molecular dioxygen or air) to be used.  
dH = Thermo(fu,5) ; % Load dH for the combustion reaction in J / mol
Tt = 298.16 ; % Load an initial temperature of 25 C
syms count ;
count = 0 ; % Start with a count of 0
syms T ; % Make T an independent variable
% Define the Cp values for the fuel, oxidant, CO2, and H2O in terms of T
% These values should end up being in J / mol K
CpFu=Thermo(fu,1)+Thermo(fu,2)*T+Thermo(fu,3)*T^2+Thermo(fu,4)*T^3;
CpOx=Thermo(ox,1)+Thermo(ox,2)*T+Thermo(ox,3)*T^2+Thermo(ox,4)*T^3;
CpCD=Thermo(6,1)+Thermo(6,2)*T+Thermo(6,3)*T^2+Thermo(6,4)*T^3;
CpWV=Thermo(7,1)+Thermo(7,2)*T+Thermo(7,3)*T^2+Thermo(7,4)*T^3;
% Next load the stoichiometric coefficients, assuming the equation
% a * fuel + b * ox = c * CO2 + d * O
a = 1 ;
if fu == 1 ; 
    b = 5 ; 
    c = 3 ; 
    d = 4 ;
end
if fu == 2 ; 
    b = 6.5 ; 
    c = 4 ; 
    d = 5 ;
end
if fu == 3 ; 
    b = 8 ; 
    c = 5 ; 
    d = 6 ;
end
if ox == 5 ; 
    b = 100/21*b ; % If using air, need larger proportion of it than O2 
end
% From this point, we will divide dHcomb into 100 equal points and, by
% evaluating the CpdT definite integrals and solving, iteratively determining the
% temperature after every 1/100 of the reaction has passed.
while count < 1 ;
    count = count + 1;
    prds = count / 100 ; % proportion of products present at the end of this iteration
    rxts = 1 - prds    ; % proportion of reactants present at the end of this iteration
    syms Tf ; % define a new independent variable Tf as the final temperature of the integration
    PrdInt = int(c*CpCD+d*CpWV,T,Tt,Tf) ; % Integral of the product temperature
    RxtInt = int(a*CpFu+b*CpOx,T,Tt,Tf) ; % Integral of the reactant temperature
    % Now evaluate how much 1 / 100 of dHcomb has raised the temperature;
    %func=prds*PrdInt+rxts*RxtInt+dH/100
    Tf = fsolve(@func, 298.16, optimset('Display','off')) ;
    Tf = solve('func = 0');
    Tt = Tf;
    % above: Set the final temperature as the initial temperature for the next iteration
end
Tk = Tt ;