clear
echo on
%Solution to Problem 21
%
%This script determines the optimum concentration of B to
%return the maximum concentration of AB.  It uses the functions
%about.m and bal.m.  The function bal.m takes in estimates of
%the concentrations of A, B, AB, and A2B and returns an array
%of zeros when the four stoichiometric and equilibrium balances
%are satisfied.  The function is called using the root finding routine
%fsolve in the function about.m to determine the equilibrium
%concentrations as a function of the initial concentration of B.  The
%function about.m returns the negative of AB, the desired product.
%Thus:

%The optimum concentration of B is:
optb=fmin('about',0,5)

%And the resulting maximum in AB is:
maxab=-about(optb)

%We can also do a little graphics:
x=[0:.1:5];
out=zeros(size(x));
for i=1:length(x)
 out(i)=about(x(i));
 echo off
end
echo on

plot(x,-out)
hold on
plot(optb,maxab,'og')
hold off
xlabel('initial concentration of B')
ylabel('equilibrium concentration of AB')


____________________________



function y=about(x)
%This function takes in a value of a0 and returns
%the negative of the concentration of ab.  It returns
%the negative because we want to use it under a function
%minimization routine to -maximize- ab.  The concentration
%of ab is determined by using a root solver on the function
%balance.m.  Note that all concentrations have to be positive,
%this can be forced by the judicious use of the absolute
%value command, or by slack variables.

global apass
apass=x;

conc=fsolve('bal',.5*[1,1,1,1]');
y=-abs(conc(3));


______________________________


function y=bal(x)
%This function takes in an array of concentrations
%and returns an array of zeros when they satisfy
%both stoichiometry and equilibrium relations.  The
%initial concentration of a is passed in through the
%global variable apass.

global apass
a0=apass;


x=abs(x);

a=x(1);
b=x(2);
ab=x(3);
a2b=x(4);

b0=1;
k1=2;
k2=5;

y=zeros(size(x));
y(1)=ab-k1*a*b;
y(2)=a2b-k2*a*ab;
y(3)=a0-(a+ab+2*a2b);
y(4)=b0-(b+ab+a2b);