clear
global k1 k2 k3 k4 k5 Cto
k1=36.96;  % in (mol-s^2/liter)^-0.5
k2=210;   % in (mol-s^2/liter)^-0.5
k3=10;   % in (mol-s^2/liter)^-0.5
k4=5.67;   % in (mol-s^2/liter)^-0.5
k5=2.1;	  % in (mol-s^2/liter)^-0.5
     % 
     %  Total Concentration set by operating 
     %  pressure (35 atm) and temp (800 K)
	 %
Cto=0.53;  % in mol/l
     %
     % Molar flow of C11 specified by Cto 
     % volumetric flow rate (1 m^3/s = 1000 liter/s)
	 %
Fo(1)=Cto*1000;  % mol/l * l/s = mol/s
     %
	 % The process engineer can dial in feed ratio of H2:C11
	 % set as 2:1 according to problem statement
	 %
Feedratio=.75;  
Fo(7) = Feedratio*Fo(1);
     %
	 % All other initial flows are zero
	 % F(1)=C11, F(2)=C10, F(3)=C9, F(4)=C8, F(5)=C7, 
	 % F(6)=C6, F(7)=H2, F(8)=CH4
	 %
Fo(2)=0;     % mol/s
Fo(3)=0;     % mol/s
Fo(4)=0;     % mol/s
Fo(5)=0;     % mol/s
Fo(6)=0;     % mol/s
Fo(8)=0;     % mol/s
Vfin=130;  % liters, pick a final volume big enough to see whole behavior
VSPAN = [0 Vfin];
[V,F]=ode45('odeP6_16',VSPAN,Fo);  %V in liters, F in mol/s
FC11=F(:,1);
FC10=F(:,2);
FC9=F(:,3);
FC8=F(:,4);
FC7=F(:,5);
FC6=F(:,6);
FH2=F(:,7);
FCH4=F(:,8);

     %  Plot the molar flow rate profiles
plot(V,FC11,'r-',V,FC10,'y-',V,FC9,'g-',V,FC8,'b-',V,FC7,'r-.',V,FC6,'y-.',V,FH2,'g-.',V,FCH4,'b-.')
xlabel('Volume (liters)')
ylabel('Molar Flow (mol/s)')
title('Problem 6-16')

     %  Plot the yield for each benzene compound
%XFC11=(Fo(1)-FC11)/Fo(1);  % This is the conversion of C11 to products
%YFC10=FC10./(Fo(1)-FC11);
%YFC9=FC9./(Fo(1)-FC11);
%YFC8=FC8./(Fo(1)-FC11);
%YFC7=FC7./(Fo(1)-FC11);
%YFC6=FC6./(Fo(1)-FC11);
%plot(V,XFC11,'w-',V,YFC10,'m-.',V,YFC9,'c--',V,YFC8,'r-',V,YFC7,'g-.',V,YFC6,'b--')
%xlabel('Volume (liters)')
%ylabel('Yields')
