## Code for Figure 10.12

### 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
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
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];
```