clear all
% set up global variables for use in the DE routine.
global q uy uz Vmax Ky Kz eps delta

t0  = 0;     %initial time
t1  = 40;    % final time
z0  = 0;     % initial ?? conc
y0  = 0;     % initial nitrogen conc
x0  = 4.9;   % initial yeast concentration

% initial state for use in DE solver.

S0  = [x0;y0;z0];

q   = 0.15; % dilution rate
eps = q/1.1;

Vmax  =  12.8909;
Ky    =   0.0499;
Kz    =   5.7694;
delta =   0.0592;
uz    = 134.6042;

uy    =  68;

[t,St] = ode45('chemostat_dynamics_2n',[t0,t1],S0);
    
subplot(311)
plot(t,St(:,1)),title('Biomass')

subplot(312)
plot(t,St(:,2)),title('nitrogen')

subplot(313)
plot(t,St(:,3)),title('mystery nutrient')