function dYdt = lotka_ode(t,Y)
% assume 'Y' is a vector of form Y = [x; y; z];

x = Y(1);
y = Y(2);
z = Y(3);

global k1 M1 k2 M2 c c2 gamma delta zeta



prey1 = x.*(k1.*(M1-x) - c.*y);
pred = -y.*(gamma - (x./delta) - zeta.*z);
prey2 = z.*(k2.*(M2 - z) - c2.*y);

dYdt = [prey1; pred; prey2];
 

     