function F = solve_ss(solution,pars)

TAU      = pars(1)   ;
BETA     = pars(2)   ;
ALPHAD0  = pars(3)   ;
ALPHAD1  = pars(4)   ;
ALPHAC0  = pars(5)   ;
ALPHAC1  = pars(6)   ;
ETAD     = pars(7)   ;
ETAC     = pars(8)   ;
PE       = pars(9)   ;   
CHI      = pars(10)  ;
ZC       = pars(11)  ;
BU       = pars(12)  ;
SIGMA    = pars(13)  ;
DELTA    = pars(14)  ;
ZD       = pars(15)  ;
DELTAH   = pars(16)  ;
GAMMA1   = pars(17)  ;
GAMMA2   = pars(18)  ;
SI       = pars(19)  ;
TE       = pars(20)  ;
S        = pars(21)  ;
p_c      = pars(22)  ;

n_d     = solution(1);
n_c     = solution(2);
u_d     = solution(3);
p_d     = solution(4);

%labor market defination
u_c = 1 - n_d - n_c - u_d;
u  = u_c + u_d;
v_d = ETAD*n_d;
v_c = ETAC*n_c;
theta_d = v_d / u_d;
theta_c = v_c / u_c;
k_d  =  ALPHAD0*theta_d^ALPHAD1;
k_c  =  ALPHAC0*theta_c^ALPHAC1;

%energy and envrionment
e = ((ZD^(SI-1)*n_d^((SI-1)/SI))/((TE+PE)^(SI-1)-ZD^(SI-1)))^(SI/(SI-1));
%h = (e^GAMMA1)/((1+1500*S)^GAMMA2);
h = e^GAMMA1;
h_row = 2.79*h;
%h_row = 4.0629; 
H = (h + h_row)/DELTAH ;
Omega = 1/(1 + H + H^2);

%Frim-energy input
y_d = ZD*Omega*(n_d^((SI-1)/SI)+e^((SI-1)/SI))^(SI/(SI-1));  
y_c = ZC*n_c;

%market clearing
c_d = y_d - k_d*v_d;
c_c = y_c - k_c*v_c;
c  = (c_d^(1-1/SIGMA) + c_c^(1-1/SIGMA))^(SIGMA/(SIGMA-1));

%price
P = (p_d^(1-SIGMA)+p_c^(1-SIGMA))^(1/(1-SIGMA));

%profit/n
pi_d = p_d*ZD*Omega*(n_d^((SI-1)/SI)+e^((SI-1)/SI))^(1/(SI-1))*n_d^(-1/SI);
pi_c = p_c*ZC;

%Wage
w_d = (1-DELTA)/(1-TAU)*(BU+c*CHI)+DELTA*(p_d*ZD*Omega*(n_d^((SI-1)/SI)+e^((SI-1)/SI))^(1/(SI-1))*n_d^(-1/SI)+p_d*theta_d*k_d);
w_c = (1-DELTA)/(1-TAU)*(BU+c*CHI)+DELTA*(p_c*ZC+p_c*theta_c*k_c);

%Firms Value
J_d = (pi_d - w_d)/(1 - BETA*(1-ETAD));
J_c = (pi_c - w_c)/(1 - BETA*(1-ETAC));

F(1) = p_d*k_d - BETA * J_d;
F(2) = p_c*k_c - BETA * J_c;
F(3) = theta_d / (1 - BETA * (1 - ETAD - theta_d))* ((1 - TAU) * w_d - CHI * c - BU)   - theta_c / (1 - BETA * (1 - ETAC - theta_c))* ((1 - TAU) * w_c - CHI * c - BU);
F(4) = P * c + S - w_d * n_d - w_c * n_c - TE * e ;

end

   
