% ---- This function is the initial guess of the ODE solution, in the form of y(x) -----
function yinit = Finitwithproduction(x)   
global rho sigma sigmasqr lambda gamma A phi r elbarnew wlbarnew qinf pinf plbarnew qlbarnew eubar

% Note: Finitwithproduction returns a linear interpolation of functions p(e) and q(e)
  pprime=(pinf-plbarnew)/(eubar-elbarnew);
  qprime=(qinf-qlbarnew)/(eubar-elbarnew);
   yinit = [ plbarnew+pprime*x
       qlbarnew+qprime*x
      pprime
      qprime];
  