function [ dSdt ] = chemostat_dynamics_2n(t,S) % this function defines the chemostat differential equation global q uy uz Vmax K M eps delta dSdt = zeros(size(S)); x = S(1); % cells y = S(2); % NH z = S(3); % Carbon ryz = (y/(K+y))*(z/(M+z)); dxdt = Vmax*ryz*x - q*x; dydt = q*(uy - y) - (1/eps)*Vmax*ryz*x; dzdt = q*(uz - z) - (1/delta)*Vmax*ryz*x; dSdt(1) = dxdt; dSdt(2) = dydt; dSdt(3) = dzdt; end