function res=musigma(x)
global solF etah kappa A delta gamma rho lambda sigma m phi
pe=deval(solF,x,1);
[prime1,prime2]=deval(solF,x,3);
peprime=prime1;
peprime2=prime2;
qe=deval(solF,x,2);
[prime1,prime2]=deval(solF,x,4);
qeprime=prime1;
qeprime2=prime2;

we=pe+qe;weprime=peprime+qeprime;weprime2=peprime2+qeprime2;
%eta=etal.*(espan<eentry)+etah.*(espan>=eentry);
i=(qe-1)/kappa;
c=A-delta-i-i^2/2*kappa;

ttheta=max(we/x,1/(1-lambda)); %checked

sigmae=(x*we*sigma*(m*ttheta-1))/(we-x*m*ttheta*weprime); %checked 

Fe=(pe+qe)/x-m*ttheta*(peprime+qeprime); %checked

Ge=c*kappa*Fe+gamma*qe*qeprime*m*(1-ttheta)*we; % checked

musigmae=c*kappa/Ge*((m*(1-ttheta)*(rho+gamma*i-gamma*(0.5*(qe*qeprime2+qeprime^2)*sigmae^2)/(c*kappa)...
    -0.5*gamma*(1+gamma)*(c*sigma-qe*qeprime*sigmae/kappa)^2/c^2)-i-etah)*we...
    +m*ttheta*(0.5*(weprime2)*sigmae^2+(A-delta)/phi+(pe-(1-phi)/phi)*i-(1-phi)/phi*kappa*0.5*(i^2)+delta*(1-qe)));  % rewritten, checked


mue=musigmae-sigma.*sigmae;

res=[mue;sigmae;i+delta];