function F=solve_steady(lp,ls,omegas,hp,pp,PI,DELTA,THETA_S,THETA_P,ALPHA,As,Ap,ETA,GAMMA,TAU_S,TAU_P,KAPPA,...
    SIGMA,MU,ZETA,RL,KM,OMEGAM,BETA_B,PHI,ALPHA_S,ALPHA_P,XI,BETA_E,BETA,EPSILON,BETA_SB,ALPHA_SB,ZETAL)
R=1/BETA*PI;
pw=(EPSILON-1)/EPSILON;

int_xi=KM/(KM-1)*(OMEGAM-OMEGAM^KM*omegas^(1-KM));
F_xi=1-(OMEGAM/omegas)^KM;
lambda_c=1-BETA_B*R/PI;

Rs=(ALPHA_S*ls+1-TAU_S-lambda_c*(1-GAMMA))/BETA_B*PI+RL*TAU_S;
Rp=(ALPHA_P*lp+1-TAU_P-lambda_c*(1-GAMMA*ZETAL))/BETA_B*PI+RL*TAU_P;
kp=lp*Rp/THETA_P/PI;
ks=ls*Rs/THETA_S/PI;
lp_bar=XI*kp;
Rp_bar=(1+lp_bar*ALPHA_SB)*PI/BETA_SB/(1-int_xi);

Xb=1-F_xi+(1-MU)*int_xi/omegas;
rb=PI/BETA_SB/(1-F_xi);
b=lp_bar*Rp_bar*omegas/rb;
%d=((1-TAU_P)*lp+(1-TAU_S)*ls+b-GAMMA*((1-TAU_P)*lp+(1-TAU_S)*ls+b*ZETA))/(1-rr);
d=lp+ls+b-GAMMA*(ls+ZETAL*lp+ZETA*b);
yp=Ap*kp^ALPHA*hp^(1-ALPHA);
wage=(1-ALPHA)*pp*yp/hp;
y=yp/(1-PHI)^SIGMA/(pp/pw)^(-SIGMA);
DIV_P=(-DELTA)*kp+pp*yp-wage*hp+(1-Rp/PI)*lp+(1-Rp_bar/PI)*lp_bar;
k=ks+kp;
i=DELTA*k;
DWL=int_xi*MU*lp_bar*Rp_bar;
ps=((pw^(1-SIGMA)-(1-PHI)^(SIGMA)*pp^(1-SIGMA))/((PHI)^SIGMA))^(1/(1-SIGMA));
ys=(ps/pw)^(-SIGMA)*PHI^(SIGMA)*y;
hs=(1-ALPHA)*ps*ys/(wage);
pi_SB=((1-int_xi)*lp_bar*Rp_bar-(1-F_xi)*rb*b)/PI-lp_bar+b-ALPHA_SB/2*lp_bar^2;
pi_CB=(Rp*lp+Rs*ls+rb*b*Xb-R*d-RL*(TAU_P*lp+TAU_S*ls))/PI-ALPHA_P/2*lp^2-ALPHA_S/2*ls^2-(1-TAU_S)*ls-(1-TAU_P)*lp-b+d;
DIV_S=(-DELTA)*ks+ps*ys-wage*hs+(1-Rs/PI)*ls;
ch=y-DWL-i;
h=(wage/ch/KAPPA)^(1/ETA);

F(1)=BETA_B*rb*Xb/PI-1+lambda_c*(1-ZETA*GAMMA);
F(2)=ys-As*ks^ALPHA*hs^(1-ALPHA);
F(3)=BETA_E*((1-DELTA)+ALPHA*ps*ys/ks-THETA_S)+THETA_S/Rs*PI-1;
F(4)=hp+hs-h;
F(5)=BETA_E*((1-DELTA)+ALPHA*pp*yp/kp-THETA_P-XI*Rp_bar/PI)+XI+THETA_P/Rp*PI-1;
