function dh=Fvecwithproduction(x,y)
global rho sigma lambda gamma A phi delta kappa etah etal eentry maversion zeta m
% Here x is e in the paper. y is a vector y(x) = [ p(x), q(x), p'(x), q'(x)]
p=y(1,:);q=y(2,:);pprime=y(3,:);qprime=y(4,:);
w = zeta * p + q;
wprime  = zeta * pprime + qprime;
theta=max(w./x,1/(1-lambda));
i=(q-1)/kappa;   % This is \hat{i}  in the notes.
c=A-delta-i-kappa.*(i.^2)./2;
F=w./x- m.*theta.*wprime;
G=c.*kappa.*F+gamma*q.*qprime.*m.*(1-theta).*w;
eta=etal.*(x<eentry)+etah.*(x>=eentry); 
sigmae=x.*w.*sigma.*(m.*theta-1)./(w-x.*m.*theta.*wprime);

a11=pprime.*c.*kappa./G.*(-m.*(1-theta).*gamma.*0.5.*q.*sigmae.^2./c./kappa.*w+m.*theta.*0.5.*sigmae.^2)...
    +p.*gamma./G.*(q.*qprime.*m.*theta.*0.5.*sigmae.^2+0.5.*F.*q.*sigmae.^2);
a12=pprime.*c.*kappa./G .*m.*theta.*0.5.*zeta.*sigmae.^2+0.5.*sigmae.^2+p.*gamma./G.*q.*qprime.*m.*theta.*0.5.*zeta.*sigmae.^2;
a21=qprime.*(c.*kappa./G.*(-m.*(1-theta).*gamma.*0.5.*q.*sigmae.^2./c./kappa.*w+0.5.*m.*theta.*sigmae.^2))...
    +0.5.*sigmae.^2+gamma.*q./G.*(q.*qprime.*m.*theta.*0.5.*sigmae.^2+0.5.*F.*q.*sigmae.^2);
a22=qprime.*c.*kappa./G.*m.*theta.*0.5.*zeta.*sigmae.^2+(q.^2).*gamma./G.*qprime.*m.*theta.*0.5.*zeta.*sigmae.^2;

mm=(m.*(1-theta).*(rho+gamma.*i-gamma.*qprime.^2.*sigmae.^2./(2.*c.*kappa)...
    -0.5.*gamma.*(1+gamma).*(c.*sigma-q.*qprime.*sigmae./kappa).^2./c.^2)-i-eta).*w ...
    + m.*theta.*((A-delta).* (zeta+phi-zeta*phi) /phi+zeta.*(p-(1-phi)/phi).*i- zeta.*(1-phi)/phi.*kappa.*0.5.*(i.^2)+delta.*(1-q));
nn=(rho+gamma.*i).*c.*kappa.*F...
    -0.5.*gamma.*F.*qprime.^2.*(sigmae.^2)-0.5.*gamma.*(1+gamma).*(F.*kappa.*(c.*sigma-q.*qprime.*sigmae./kappa).^2)./c ...
    -gamma.*q.*qprime.*((-i-eta).*w+m.*theta.*((A-delta).*(zeta+phi-zeta*phi) /phi+zeta.*(p-(1-phi)/phi).*i-zeta.*(1-phi)/phi.*kappa.*0.5.*(i.^2)+delta.*(1-q))) ;

% A big change:  (zeta+phi-zeta*phi) /phi, instead of 1/phi

b1=(p.*sigma+pprime.*sigmae).*sigma.*maversion.*theta.*((w-x.*wprime)./(x.*F))...
    -pprime.*c.*kappa./G.*(mm)...
    -(1-phi)/phi.*c-p.*i+p./G.*nn;
b2=(q.*sigma+qprime.*sigmae).*sigma.*maversion.*theta.*((w-x.*wprime)./(x.*F))...
    -qprime.*c.*kappa./G.*(mm)...
    -A+q.*delta+q./G.*nn;

det=a11.*a22-a12.*a21;

q2prime=(a22.*b1-a12.*b2)./det;
p2prime=(-a21.*b1+a11.*b2)./det;

dh=[y(3,:)
    y(4,:)
    p2prime
    q2prime];