syms w_t2 w_t1 z_t2 z_t1 py_t1 py_t2 theta K

f = mu * theta^(1 - eta);
q = f / theta;

uuuM=(F_t2 * d) / (d + f);
eeeM= F_t2 - uuuM;

uuuN=(F_t1 * d)/(d + h * f);
eeeN= F_t1 - uuuN;

uuu=uuuM + uuuN;

Aj  = rr + bbeta * d;
Au_t2 = Aj + bbeta * f;
Au_t1 = Aj + bbeta * h * f;

J_t2 = (py_t2 - w_t2) / Aj;
J_t1 = (py_t1 - w_t1) / Aj;
W_U_t2 = (w_t2 * (1 - t) - z_t2) / Au_t2;
W_U_t1 = (w_t1 * (1 - t) - z_t1) / Au_t1;

L   = Z * ((chi * ((eeeN)^((epsilon - 1) / epsilon)) + (1 - chi) * ((eeeM)^((epsilon - 1) / epsilon)))^(epsilon / (epsilon - 1)));
  Y   = (K^aalpha) * (L^(1 - aalpha)); % F(L,K)

%% Steady State

eq1_t2_SS = (1 - t) * b * J_t2 == (1 - b) * W_U_t2; 
eq1_t1_SS = (1 - t) * b * J_t1 == (1 - b) * W_U_t1; 

eq2_t1_SS = py_t2  == (Z^(1 - 1 / epsilon)) * (1 - aalpha) * (Y / L) * (1 - chi) * ((L / eeeM)^(1 / epsilon));
eq2_t2_SS = py_t1  == (Z^(1 - 1 / epsilon)) * (1 - aalpha) * (Y / L) * chi * ((L / eeeN)^(1 / epsilon));

eq3_t2_SS = z_t2 == 0.612 * w_t2 * (1 - t);
eq3_t1_SS = z_t1 == 0.612 * w_t1 * (1 - t);

eq4_SS= c == bbeta * q * ((uuuM / uuu) * J_t2...
                   + h *  (uuuN / uuu) * J_t1);

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

eqs_SS=[eq1_t2_SS, eq1_t1_SS, eq2_t2_SS, eq2_t1_SS, ...
        eq3_t2_SS, eq3_t1_SS, eq4_SS, eq5_SS];

guess_SS=xguess;

vars_SS=[w_t2, w_t1, z_t2, z_t1, py_t1, py_t2, theta, K];
 
        [w_t2, w_t1, z_t2, z_t1, py_t1, py_t2, theta, K]=vpasolve(eqs_SS, vars_SS, guess_SS);

xguess=[w_t2, w_t1, z_t2, z_t1, py_t1, py_t2, theta, K];

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