function res=Fbcwithproduction(ya,yb);
global rho sigma sigmasqr lambda gamma A phi r plbarnew qlbarnew pinf qinf
% ya = [  p(e_), q(e_), p'(e_), q'(e_)  ];   yb = [  p(e^-), q(e^-), p'(e^-), q'(e^-)  ]
%B=phi*(R/r+cstar)*(theta+cstar)/(r*theta+phi*theta+phi*cstar-R);
% res=[ya(1)-plbarnew
%     ya(2)-qlbarnew
%     yb(1)-pinf
%     yb(2)-qinf];
% res=[ya(2)
%     yb(1)-winf];
res=[ya(1)-plbarnew
    ya(2)-qlbarnew
    yb(3)-0
    yb(4)-0];
