function ff = solvess(habit,betap,betaip,ksih,LTV,alpha,deltak,deltah,iota,betab,gammab,rho_b)

rd=1/betap-1;
rk=1/betap-1+deltak;
rl=1/betab-1;
qk=1;
qh=1;
N1=((1-betaip*(1-deltah))-LTV*((1-gammab)/betab+gammab/betap-1/betaip)*(betaip/(betaip*rho_b-1)));

eps1=ksih*(1-habit)/((1-betap*(1-deltah))*(1-habit*betap));
eps2=ksih*(1-habit)/((1-habit*betaip)*N1);

syms qhyh_y;
ci_y=(1-alpha)*iota/(1+(deltah+rl*(1-habit)*ksih/((1-habit*betaip)*N1))*eps2);
cb_y=(rl-rd*gammab)*LTV/(1-rho_b)*eps2*ci_y;
ik_y=deltak*alpha/rk;
ih_y=deltah*qhyh_y;
cp_y=1-ci_y-cb_y-ik_y-ih_y;

qhhp_y=eps1*cp_y;
qhhi_y=eps2*ci_y;

eq=(qhhp_y+qhhi_y)-qhyh_y;
qhyh_y=solve(eq,qhyh_y);
ff=double(qhyh_y);