Figure 10.12 (page 598):

The two steady-state biomass and substrate concentrations versus dilution rate; stable (solid), unstable (dashed).

Code for Figure 10.12

Text of the GNU GPL.

main.m

% straighten out the chemostat
%

global mum D Ks Sf y
mum = 1;
D = 1.05;
Ks = 1;
Sf = 5;
y = 1;


Dopt = mum*(1-sqrt(Ks/(Ks+Sf)))
Dc = mum*Sf/(Ks+Sf)

nDs = 25;
Dvec1 = linspace (0.001, Dc, nDs)';
for i = 1:nDs
  D = Dvec1(i);
  %compute eigenvalues of jacobian for both steady states
  %first steady state is washout
  xs11(i,:) = [0, Sf];
  Ss = D*Ks/(-D+mum);
  xs21(i,:) = [y*(Sf-Ss), Ss];
  lam1(i,:) = sort(eig(jacob(xs11(i,:))));
  lam2(i,:) = sort(eig(jacob(xs21(i,:))));
end


Dmax = Dc;
Dvec2 = linspace (Dc, 0.99, nDs)';
for i = 1:nDs
  D = Dvec2(i);
  %compute eigenvalues of jacobian for both steady states
  %first steady state is washout
  xs12(i,:) = [0, Sf];
  Ss = D*Ks/(-D+mum);
  xs22(i,:) = [y*(Sf-Ss), Ss];
  lam12(i,:) = sort(eig(jacob(xs12(i,:))));
  lam22(i,:) = sort(eig(jacob(xs22(i,:))));
end


% compute production rate
DX1 = Dvec1.*xs21(:,1);
DX2 = Dvec2.*xs22(:,1);

figure(1)
plot(Dvec1, [xs11(:,1), xs21(:,1), DX1], ...
     Dvec2, [xs12(:,1), xs22(:,1), DX2], '-')
axis([0,1,-6,6])

figure(2)
plot(Dvec1, [xs11(:,2), xs21(:,2)], ...
     Dvec2, [xs12(:,2), xs22(:,2)], '-')
axis([0,1,-0.5,7])

table1 = [Dvec1, xs11, xs21, DX1];
table2 = [Dvec2, xs12, xs22, DX2];

save "chemoss.dat"  table1 table2

jacob.m

function J = jacob(x)
  global mum D Ks Sf y
  X = x(1);
  S = x(2);
  mug = (mum*S)/(Ks+S);
  tmp = Ks*mum/(Ks+S)^2;
  J = [mug - D, X*tmp; -mug/y, -D-X/y*tmp];