syms w_L_t2 w_H_t2 w_L_t1 w_H_t1...
     z_L_t2 z_H_t2 z_L_t1 z_H_t1...
     py_L_t2 py_L_t1 py_H_t1 py_H_t2...
     theta_L theta_H K

f_H = mu_H * theta_H^(1 - eta);
f_L = mu_L * theta_L^(1 - eta);
q_L = f_L / theta_L; 
q_H = f_H / theta_H;

uuuH_t2=(F_H_t2 * d_H) / (d_H + f_H);
eeeH_t2= F_H_t2 - uuuH_t2;
uuuL_t2=(F_L_t2 * d_L) / (d_L + f_L);
eeeL_t2= F_L_t2 - uuuL_t2;

uuuH_t1=(F_H_t1 * d_H)/(d_H + hH_t1 * f_H);
eeeH_t1= F_H_t1 - uuuH_t1;
uuuL_t1=(F_L_t1 * d_L)/(d_L + hL_t1 * f_L);
eeeL_t1= F_L_t1 - uuuL_t1;

uuuL=uuuL_t2 + uuuL_t1;
uuuH=uuuH_t2 + uuuH_t1;

Aj_L  = rr + bbeta * d_L;
Aj_H  = rr + bbeta * d_H;
Au_H_t2 = Aj_H + bbeta * f_H;
Au_H_t1 = Aj_H + bbeta * hH_t1 * f_H;
Au_L_t2 = Aj_L + bbeta * f_L;
Au_L_t1 = Aj_L + bbeta * hL_t1 * f_L;

J_H_t2 = (py_H_t2 - w_H_t2) / Aj_H;
J_H_t1 = (py_H_t1 - w_H_t1) / Aj_H;
J_L_t2 = (py_L_t2 - w_L_t2) / Aj_L;
J_L_t1 = (py_L_t1 - w_L_t1) / Aj_L;
W_U_H_t2 = (w_H_t2 * (1 - t) - z_H_t2) / Au_H_t2;
W_U_H_t1 = (w_H_t1 * (1 - t) - z_H_t1) / Au_H_t1;
W_U_L_t2 = (w_L_t2 * (1 - t) - z_L_t2) / Au_L_t2;
W_U_L_t1 = (w_L_t1 * (1 - t) - z_L_t1) / Au_L_t1;

l_L   = Z_L * (eeeL_t2 + eeeL_t1 * lambda_L);
l_H   = Z_H * ((chi * ((eeeH_t1)^((epsilon - 1) / epsilon)) + (1 - chi) * ((eeeH_t2)^((epsilon - 1) / epsilon)))^(epsilon / (epsilon - 1)));
  L   = (phi * ((l_L)^((rho - 1) / rho)) + (1 - phi) * ((l_H)^((rho - 1) / rho)))^(rho / (rho - 1));
  Y   = (K^aalpha) * (L^(1 - aalpha)); % F(L,K)

%% Steady State

eq1_H_t2_SS = (1 - t) * b_H * J_H_t2 == (1 - b_H) * W_U_H_t2; 
eq1_H_t1_SS = (1 - t) * b_H * J_H_t1 == (1 - b_H) * W_U_H_t1; 
eq1_L_t2_SS = (1 - t) * b_L * J_L_t2 == (1 - b_L) * W_U_L_t2; 
eq1_L_t1_SS = (1 - t) * b_L * J_L_t1 == (1 - b_L) * W_U_L_t1; 

eq2_L_t1_SS = py_L_t2  == Z_L * (1 - aalpha) * (Y / L) * phi * ((L / l_L)^(1 / rho));
eq2_H_t1_SS = py_H_t2  == (Z_H^(1 - 1 / epsilon)) * (1 - aalpha) * (Y / L) * (1 - phi) * (1 - chi) * ((L / l_H)^(1 / rho)) * ((l_H / eeeH_t2)^(1 / epsilon));
eq2_L_t2_SS = py_L_t1  == lambda_L * py_L_t2;
eq2_H_t2_SS = py_H_t1  == (Z_H^(1 - 1 / epsilon)) * (1 - aalpha) * (Y / L) * (1 - phi) * chi * ((L / l_H)^(1 / rho)) * ((l_H / eeeH_t1)^(1 / epsilon));

eq3_H_t2_SS = z_H_t2 == 0.612 * w_H_t2 * (1 - t);
eq3_H_t1_SS = z_L_t2 == 0.612 * w_L_t2 * (1 - t);
eq3_L_t2_SS = z_H_t1 == 0.612 * w_H_t1 * (1 - t);
eq3_L_t1_SS = z_L_t1 == 0.612 * w_L_t1 * (1 - t);

eq4_H_SS= cH == bbeta * q_H * ((uuuH_t2 / uuuH) * J_H_t2...
                      + hH_t1 *  (uuuH_t1 / uuuH) * J_H_t1);

eq4_L_SS= cL == bbeta * q_L * ((uuuL_t2 / uuuL) * J_L_t2...
                      + hL_t1 *  (uuuL_t1 / uuuL) * J_L_t1);

eq5_SS= aalpha * (Y / K) == r + ddd;

eqs_SS=[eq1_H_t2_SS, eq1_H_t1_SS, eq1_L_t2_SS, eq1_L_t1_SS, ...
        eq2_H_t2_SS, eq2_H_t1_SS, eq2_L_t2_SS, eq2_L_t1_SS, ...
        eq3_H_t2_SS, eq3_H_t1_SS, eq3_L_t2_SS, eq3_L_t1_SS, ...
        eq4_H_SS, eq4_L_SS, eq5_SS];

guess_SS=xguess;

vars_SS=[w_L_t2, w_H_t2, w_L_t1, w_H_t1...
         z_L_t2, z_H_t2, z_L_t1, z_H_t1...
         py_L_t2, py_L_t1, py_H_t1, py_H_t2...
         theta_L, theta_H, K];
 
        [w_L_t2, w_H_t2, w_L_t1, w_H_t1,...
         z_L_t2, z_H_t2, z_L_t1, z_H_t1,...
         py_L_t2, py_L_t1, py_H_t1, py_H_t2,...
         theta_L, theta_H, K]=vpasolve(eqs_SS, vars_SS, guess_SS);

xguess=[w_L_t2, w_H_t2, w_L_t1, w_H_t1, ... 
        py_L_t2, py_L_t1, py_H_t1, py_H_t2, ...
        z_L_t2, z_H_t2, z_L_t1, z_H_t1, ...
        theta_H, theta_L, K];

save('stst_guess.mat','xguess')