function dydt = odefun12(t,w)
x = w(1);
y = w(2);
z1 = w(3);
z2 = w(4);

qin=1.5;c=1;

if isreal(z2)==1 && isreal(z1)==1
    s=0; % 3 real roots of f3
else 
    s=1; % 1 real root
end  
% if Q<2.6919
%     s=0;
% else
%     s=1;
% end
f1 = qin -c*x-s*z1-(1-s)*z2;      % differential eq
f2 = 1/y;                         % differential eq 
f3 = (z1^3 + y*z1^2+z1-1.2);      % alg eq
f4 = (z2^3 + y*z2^2+z2-1.2);      % alg eq, can also be written with s*

dydt = [f1;f2;f3;f4];
end