function [ dSdt ] = chemostat_dynamics_2n(t,S)
% this function defines the chemostat differential equation
global q uy uz V Ky Kz eps delta

dSdt = zeros(size(S));

x = S(1); % cells
y = S(2); % Nut1
z = S(3); % Nut2

ryz = (y/(Ky+y))*(z/(Kz+z));

dxdt  = V*ryz*x - q*x;
dydt  = q*(uy - y) - (1/eps)*V*ryz*x;
dzdt  = q*(uz - z) - (1/delta)*V*ryz*x;

dSdt(1) = dxdt;
dSdt(2) = dydt;
dSdt(3) = dzdt;
end