function dx = chemostatdynamics(t,x)
global q u V k r 

dx = zeros(size(x));

dx(1) = (q*u) - (q*x(1)) - (x(2)*((V*x(1))/(x(1)+k)));
dx(2) = (-q*x(2)) + ((r*x(2))*((V*x(1))/(x(1)+k)));

end