function dx = nutrientPop(t,x)
global a b epsiLon g

dx = zeros(size(x));

dx(1) = a*x(1)*x(2) - b*x(1);
dx(2) = -a*epsiLon*x(1)*x(2) + g*x(2);

end