syms      pyH_1 pyL_1 pyH_2 pyL_2 wL_1 wL_2 wH_1 wH_2...
          theta_H theta_L muL_1 muH_1 muL_2 muH_2 p2...
          XL_1 XL_2 XH_1 XH_2 JL_1 JL_2 JH_1 JH_2

%% Auxiliares

mL_2=theta_L^(1-eta)*muL_2;
mH_2=theta_H^(1-eta)*muH_2;
mL_1=theta_L^(1-eta)*muL_1;
mH_1=theta_H^(1-eta)*muH_1;
mL=mL_2+mL_1;
mH=mH_2+mH_1;
q_L=mL/theta_L;
q_H=mH/theta_H;

uuuH=(d*ff_H)/(d+mH);
eeeH_1=(mH_1*uuuH)/d;
eeeH_2=ff_H-uuuH-eeeH_1;

uuuL=(d*ff_L)/(d+mL);
eeeL_1=(mL_1*uuuL)/d;
eeeL_2=ff_L-uuuL-eeeL_1;

Y_2=(alpha_2*((eeeL_2)^((gggamma-1)/gggamma))+(1-alpha_2)*((eeeH_2)^((gggamma-1)/gggamma)))^(gggamma/(gggamma-1));
Y_1=(alpha_1*((eeeL_1)^((gggamma-1)/gggamma))+(1- alpha_1)*((eeeH_1)^((gggamma-1)/gggamma)))^(gggamma/(gggamma-1));

Y= p2*Y_2/(Y_1*p1+Y_2*p2);
P=(p1*(p1/a)^(-sssigma)+p2*(p2/(1-a))^(-sssigma))^(1/(1-sssigma));

% Present value of Workers
eq1=XL_1==((wL_1*(1-t)-z_L)/P-(mL_2*XL_2))/(r+d+mL_1);
eq2=XL_2==((wL_2*(1-t)-z_L)/P-(mL_1*XL_1))/(r+d+mL_2);
eq3=XH_1==((wH_1*(1-t)-z_H)/P-(mH_2*XH_2))/(r+d+mH_1);
eq4=XH_2==((wH_2*(1-t)-z_H)/P-(mH_1*XH_1))/(r+d+mH_2);

% Present value of Firms
eq5=JL_1==(pyL_1-wL_1+d*VL_1)/(d+r);
eq6=JH_1==(pyH_1-wH_1+d*VH_1)/(d+r);
eq7=JL_2==(pyL_2-wL_2+d*VL_2)/(d+r);
eq8=JH_2==(pyH_2-wH_2+d*VH_2)/(d+r);

eq9= cL_1+(r+q_L)*VL_1==q_L*JL_1;
eq10=cL_2+(r+q_L)*VL_2==q_L*JL_2;
eq11=cH_1+(r+q_H)*VH_1==q_H*JH_1;
eq12=cH_2+(r+q_H)*VH_2==q_H*JH_2;

% Nash Barganing equations
eq13=XL_1*P*(1-bL_1)==(JL_1-VL_1)*(1-t)*bL_1;
eq14=XL_2*P*(1-bL_2)==(JL_2-VL_2)*(1-t)*bL_2;
eq15=XH_1*P*(1-bH_1)==(JH_1-VH_1)*(1-t)*bH_1;
eq16=XH_2*P*(1-bH_2)==(JH_2-VH_2)*(1-t)*bH_2;

% FOC of each worker type in each sector
eq17=pyL_1==p1*alpha_1*((Y_1/eeeL_1)^(1/gggamma));
eq18=pyH_1==p1*(1-alpha_1)*((Y_1/eeeH_1)^(1/gggamma));
eq19=pyL_2==p2*alpha_2*((Y_2/eeeL_2)^(1/gggamma));
eq20=pyH_2==p2*(1-alpha_2)*((Y_2/eeeH_2)^(1/gggamma));

eq21=Y==p2*(p2/(1-a))^(-sssigma)/(p1*(p1/a)^(-sssigma)+p2*(p2/(1-a))^(-sssigma));
eq22=0==1-muL_2-muL_1;
eq23=0==1-muH_2-muH_1;

eqs=[eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,...
     eq12,eq13,eq14,eq15,eq16,eq17,eq18,eq19,eq20,...
     eq21,eq22,eq23];

vars=[theta_L,theta_H,muL_1,muH_1,muL_2,muH_2,p2,...
      pyL_2,wL_2,pyL_1,wL_1,pyH_1,wH_1,pyH_2,wH_2,...
      XL_1,XL_2,XH_1,XH_2,JL_1,JL_2,JH_1,JH_2];
   
range=xguess;

     [theta_L,theta_H,muL_1,muH_1,muL_2,muH_2,p2,...
      pyL_2,wL_2,pyL_1,wL_1,pyH_1,wH_1,pyH_2,wH_2,...
      XL_1,XL_2,XH_1,XH_2,JL_1,JL_2,JH_1,JH_2]=vpasolve(eqs,vars,range);