function [ dSdt ] = chemostat_2n(t,S)
% this function defines the chemostat differential equation
global q u r K V u2
dSdt = zeros(size(S));

c1 = S(1);
y = S(2);
c2 = S(3);

dc1dt = q*u - q*c1 -((y*c1*V)/(K+c1))*(c2/(c2+K));
dc2dt = q*u2 - q*c2 -((y*c1*V)/(K+c1))*(c2/(c2+K));

dydt = (y*r)*(V*c1)/(K+c1)*(c2/(c2+K)) - q*y;

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