clear all
t0 = 0;   %initial time
t1 = 1;   % final time
c0 = 0;  % initial nutrient conc
N0 = 30;
S0 = [c0;N0];
%
% the global statement below allows
%the chemostat DE function to know
%all the rate parameters
%
global q u r K V
q = 0.15
u = 30
r = 1.25
K = 5
V = 0.5

cstar = K*q/(r*V-q)
nstar = (q*u-q*cstar)/(V*cstar/(K+cstar))

[t,St] = ode45('chemostat',[t0,t1],S0);
plot(t,St)
legend('nutrient','cells')