function [ dSdt ] = chemostat_dynamics(t,S)
% this function defines the chemostat differential equation
global q un uc Kn Kc Vn Vc r

dSdt = zeros(size(S));

c1 = S(1); % c1 = carbon
c2 = S(2); % c2 = nitrogen
y  = S(3); % c3 = yeast

dc1dt = q*uc- q*c1 -((y*c1*Vc)/(Kc+c1))*(c2/(c2+Kn));
dc2dt = q*un - q*c2 -((y*c1*Vn)/(Kc+c1))*(c2/(c2+Kn));

dydt  = (y*r)*(c1)/(Kc+c1)*(c2/(c2+Kn)) - q*y;

dSdt(1) = dc1dt;
dSdt(2) = dc2dt;
dSdt(3) = dydt;
end