function ff = solvess(BETAp,BETAi,BETAr,BETAb,XI,ZETAh,TAU,DELTAh,LTV,DELTAk,RHOb,CAR,ALPHA,PSIi,PSIr,g_ss,epsilon)

rh=1+TAU-BETAp*(1-DELTAh);
rk=(1-BETAp*(1-DELTAk))/BETAp;
rd=1/BETAp-1;
qh=1;
qk=1;
rl=CAR/BETAb+(1-CAR)/BETAp-1;

eps1=ZETAh/(1+TAU-BETAi*(1-DELTAh)+LTV/(1-BETAi*RHOb)*(CAR*BETAi/BETAb+(1-CAR)*BETAi/BETAp-1));
eps2=eps1*(1-XI)/(1-XI*BETAi);

mc=(epsilon-1)/epsilon;

syms qhyh_y;
ci_y=(1-ALPHA)*PSIi*mc/(1+DELTAh*eps2+rl*LTV*(1-XI)*eps1/((1-RHOb)*(1-XI*BETAi)));
cr_y=((1-ALPHA)*PSIr*mc)/(1+ZETAh*(1-XI)/(1-XI*BETAr));
cb_y=(rl-rd*(1-CAR))*LTV*(1-XI)*eps1/((1-RHOb)*(1-XI*BETAi))*ci_y;
ik_y=DELTAk*ALPHA*mc/rk;
ih_y=DELTAh*qhyh_y;
cp_y=1-ci_y-cr_y-cb_y-ik_y-ih_y-g_ss;

qhhp_y=(1-XI)*ZETAh/((1-BETAp*(1-DELTAh))*(1-XI*BETAp))*cp_y;
qhhi_y=eps2*ci_y;
qhhr_y=ZETAh*(1-XI)/((1-XI*BETAr)*rh)*cr_y;

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