function [P_A,P_B,W_A,W_B,C_A,C_B,L_A,L_B] = Multi_ss(GAMMA_AA, GAMMA_BB, GAMMA_BA, GAMMA_AB, ALPHA_A, ALPHA_B, BETA_A, BETA_B, OMEGA, MU, DELT, DELTT, PHI_A, PHI_B)

f=@(ss)([
	%C_A = (P_A/P)^(-MU)*PHI_A*C;
    - ss(5) + (ss(1)/((PHI_A*ss(1)^(-(1-MU))+ PHI_B*ss(2)^(-(1-MU)))^(-1/(1-MU))))^(-MU)*PHI_A*(PHI_A^(1/MU)*ss(5)^((MU-1)/(MU))+ PHI_B^(1/MU)*ss(6)^((MU-1)/(MU)))^((MU)/(MU-1));
	%C_B = (P_B/P)^(-MU)*PHI_B*C;
	- ss(6) + (ss(2)/((PHI_A*ss(1)^(-(1-MU))+ PHI_B*ss(2)^(-(1-MU)))^(-1/(1-MU))))^(-MU)*PHI_B*(PHI_A^(1/MU)*ss(5)^((MU-1)/(MU))+ PHI_B^(1/MU)*ss(6)^((MU-1)/(MU)))^((MU)/(MU-1));
   
	%L_A = (W_A/W)^(1/OMEGA)*L; 
	- ss(7)+ (ss(3)/((ss(3)^((1+OMEGA)/OMEGA)+ ss(4)^((1+OMEGA)/OMEGA))^(OMEGA/(1+OMEGA))))^(1/OMEGA)*(ss(7)^(1+OMEGA)+ ss(8)^(1+OMEGA))^(1/(1+OMEGA)) ; 
 	%L_B = (W_B/W)^(1/OMEGA)*L;
	- ss(8)+ (ss(4)/((ss(3)^((1+OMEGA)/OMEGA)+ ss(4)^((1+OMEGA)/OMEGA))^(OMEGA/(1+OMEGA))))^(1/OMEGA)*(ss(7)^(1+OMEGA)+ ss(8)^(1+OMEGA))^(1/(1+OMEGA));
    
	% 1 = (DELT*ALPHA_A/(1-DELT*(1-DELTT)))^ALPHA_A*(BETA_A*ss(1)/ss(3))^(BETA_A)*GAMMA_AA^GAMMA_AA*(GAMMA_BA*ss(1)/ss(2))^GAMMA_BA);
	- 1 + (DELT*ALPHA_A/(1-DELT*(1-DELTT)))^ALPHA_A*(BETA_A*ss(1)/ss(3))^(BETA_A)*GAMMA_AA^GAMMA_AA*(GAMMA_BA*ss(1)/ss(2))^GAMMA_BA;
	% 1 = (DELT*ALPHA_A/(1-DELT*(1-DELTT)))^ALPHA_A*(BETA_A*ss(1)/ss(3))^(BETA_A)*GAMMA_AA^GAMMA_AA*(GAMMA_BA*ss(1)/ss(2))^GAMMA_BA);
	- 1 + (DELT*ALPHA_B/(1-DELT*(1-DELTT)))^ALPHA_B*(BETA_B*ss(2)/ss(4))^(BETA_B)*GAMMA_BB^GAMMA_BB*(GAMMA_AB*ss(2)/ss(1))^GAMMA_AB;
    
    % W_A*L_A/(BETA_A*P_A) = C_A + DELTT*(DELT*ALPHA_A/(1-DELT*(1-DELTT)))*W_A*L_A/(BETA_A*P_A)+GAMMA_AA*W_A*L_A/(BETA_A*P_A)+GAMMA_AB*P_B/P_A*W_B*L_B/(BETA_B*P_B);
    - ss(3)*ss(7)/(BETA_A*ss(1)) + ss(5) + DELTT*(DELT*ALPHA_A/(1-DELT*(1-DELTT)))*ss(3)*ss(7)/(BETA_A*ss(1))+GAMMA_AA*ss(3)*ss(7)/(BETA_A*ss(1))+GAMMA_AB*ss(2)/ss(1)*ss(4)*ss(8)/(BETA_B*ss(2));
    % W_A*L_A/(BETA_A*P_A) = C_A + DELTT*(DELT*ALPHA_A/(1-DELT*(1-DELTT)))*W_A*L_A/(BETA_A*P_A)+GAMMA_AA*W_A*L_A/(BETA_A*P_A)+GAMMA_AB*P_B/P_A*W_B*L_B/(BETA_B*P_B);
    - ss(4)*ss(8)/(BETA_B*ss(2)) + ss(6) + DELTT*(DELT*ALPHA_B/(1-DELT*(1-DELTT)))*ss(4)*ss(8)/(BETA_B*ss(2))+GAMMA_BB*ss(4)*ss(8)/(BETA_B*ss(2))+GAMMA_BA*ss(1)/ss(2)*ss(3)*ss(7)/(BETA_A*ss(1));
    ]);
x0=[1,1,1,1,1,1,1,1];
options=optimset('display','off','MaxFunEvals',3000,'MaxIter',3000);
rez=fsolve(f,x0,options);


P_A = rez(1);
P_B = rez(2);
W_A = rez(3);
W_B = rez(4);
C_A = rez(5);
C_B = rez(6);
L_A = rez(7);
L_B = rez(8);

end