clear all
% set up global variables for use in the DE routine.
global q uy rM K d eps

t0  = 0;     %initial time
t1  = 100;    % final time
x0  = 10;     % initial cell population
y0  = 12;    % initial yeast concentration
%
% initial state for use in DE solver.
%
S0  = [x0;y0];

%
% rate constants and problem parameters
%
rM = 5;
K  = 15;
d  =  1.8;
eps = 0.1;

uy  = 15;
q   = 0.15;

%
% call the matlab function to solve the DE for us.
%
[t,St] = ode45('chemostat_dynamics',[t0 t1],S0);
%
% Plot the results.
%
plot(t,St,'LineWidth',3)

legend('cells','food')
xlabel('time')
ylabel('states')