function F = tese_steadystate_helper(x, beta_p, beta_i, beta_e, piss, m_ess, m_ihss, m_iwss, deltak, deltah, deltakb, jh, eps_lss, phi, eps_yss, mu, alphal, alphak, eps_d, eps_biw, eps_bih, eps_be, phi_dih, phi_diw, phi_de)

r_d = piss/beta_p - 1;
r_pol = ((eps_d - 1)/eps_d)*r_d;
r_biw = (eps_biw/(eps_biw - 1))*r_pol;
r_bih = (eps_bih/(eps_bih - 1))*r_pol;
r_be = (eps_be/(eps_be - 1))*r_pol;
x_ss = eps_yss/(eps_yss - 1);
qh = 1;
qk = 1;

f1 = - x(1) + 1/x(7);
f2 = - x(2) + 1/x(8);
f3 = - x(3) + 1/x(9);
f4 = - x(4) + (1/(1 + r_bih) - beta_i/piss)*x(2);
f5 = - x(5) + (1/(1 + r_biw) - beta_i/piss)*x(2);
f6 = - x(6) + (1/(1 + r_be) - beta_e/piss)*x(3);
f7 = - x(7) + x(16)*x(21) + ((1 + r_d)/piss - 1)*(x(13) + x(14) + x(15) - x(23)) + (1 - 1/x_ss)*x(19) - deltah*x(10);
f8 = - x(8) + x(17)*x(22) - ((1 + r_bih)/piss - 1)*(phi_dih^-1)*x(13) - ((1 + r_biw)/piss - 1)*(phi_diw^-1)*x(14) - deltah*x(11);
f9 = - x(9) + x(18)/x_ss - x(16)*x(21) - x(17)*x(22) - ((1 + r_be)/piss - 1)*(phi_de^-1)*x(15) - deltak*x(20) - deltah*x(12);
f10 = - x(10) + jh/((1 - (1 - deltah)*beta_p)*x(1));
f11 = - x(11) + jh/((1 - (1 - deltah)*beta_i)*x(2) - (1 - deltah)*x(4)*m_ihss*piss);
f12 = - x(12) + jh/((1 - (1 - deltah)*beta_e)*x(3) - (1 - deltah)*x(6)*m_ess*piss);
f13 = - x(13) + (1 / (1 + r_bih))*m_ihss*piss*qh*(1 - deltah)*x(11) - (phi_dih^-1)*x(13);
f14 = - x(14) + (1 / (1 + r_biw))*m_iwss*piss*x(17)*x(22) - (phi_diw^-1)*x(14);
f15 = - x(15) + (1 / (1 + r_be))*m_ess*piss*(qk*(1 - deltak)*x(20) - qh*(1 - deltah)*x(12)) - (phi_de^-1)*x(15);
f16 = - x(16) + (eps_lss/(eps_lss - 1))*(x(21)^phi)/x(1);
f17 = - x(17) + (eps_lss/(eps_lss - 1))*(x(22)^phi)/(x(2) + x(5)*m_iwss);
f18 = - x(18) + (x(20)^alphak)*(x(21)^(mu*alphal))*(x(22)^((1 - mu)*alphal));
f19 = - x(19) + x(7) + x(8) + x(9) + deltah*(x(10) + x(11) + x(12)) + deltak*x(20);
f20 = - x(3)*qk + x(6)*m_ess*qk*piss*(1 - deltak) + beta_e*x(3)*(alphak*x(18)/(x_ss*x(20)) + qk*(1 - deltak));
f21 = - x(21) + (mu*alphal/x_ss)*(x(18)/x(16));
f22 = - x(22) + ((1 - mu)*alphal/x_ss)*(x(18)/x(17));
f23 = - x(23) + (1/(piss - 1 + deltakb))*(r_bih*x(13) + r_biw*x(14) + r_be*x(15) - r_d*(x(13) + x(14) + x(15) - x(23)));

F = [f1; f2; f3; f4; f5; f6; f7; f8; f9; f10; f11; f12; f13; f14; f15; f16; f17; f18; f19; f20; f21; f22; f23];

x0 = [1; 1;	1;	1;	1; 1; 1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1;	1];

options=optimset('MaxFunEvals',50000000,'MaxIter',50000000);

[x,fval] = fsolve(@tese_steadystate_helper,x0,options);

r_d = piss/beta_p - 1;
r_pol = ((eps_d - 1)/eps_d)*r_d;
r_biw = (eps_biw/(eps_biw - 1))*r_pol;
r_bih = (eps_bih/(eps_bih - 1))*r_pol;
r_be = (eps_be/(eps_be - 1))*r_pol;

% Dicionário:
% x(1): lam_p
% x(2): lam_i
% x(3): lam_e
% x(4): s_ih
% x(5): s_iw
% x(6): s_e
% x(7): c_p
% x(8): c_i
% x(9): c_e
% x(10): h_p
% x(11): h_i
% x(12): h_e
% x(13): b_ih
% x(14): b_iw
% x(15): b_e
% x(16): w_p
% x(17): w_i
% x(18): y_e
% x(19): y
% x(20): k
% x(21): l_p
% x(22): l_i
% x(23): k_b

disp(' ')
PIBo = x(19);
disp(['Produto otica da oferta = ' num2str(PIBo)])
disp(' ')
PIBd = x(9) + x(8) + x(7) + deltah*(x(10) + x(11) + x(12)) + deltak*x(20);
disp(['Produto otica da demanda = ' num2str(PIBd)])
disp(' ')
C = x(9) + x(8) + x(7);
disp(['Consumo agregado = ' num2str(C)])
disp(' ')
C_PIBd = C/PIBd;
disp(['Consumo/PIBd = ' num2str(C_PIBd)])
disp(' ')
I = deltak*x(20) + deltah*(x(10) + x(11) + x(12));
disp(['Investim_essnto agregado = ' num2str(I)])
disp(' ')
I_PIBd = I/PIBd;
disp(['Investim_essnto/PIBd = ' num2str(I_PIBd)])
disp(' ')
Ik = deltak*x(20);
disp(['Investim_essnto em capital = ' num2str(Ik)])
disp(' ')
Ik_PIBd = Ik/PIBd;
disp(['Investim_essnto em capital/PIBd = ' num2str(Ik_PIBd)])
disp(' ')
Ih = deltah*(x(10) + x(11) + x(12));
disp(['Investim_essnto em Imoveis = ' num2str(Ih)])
disp(' ')
Ih_PIBd = Ih/PIBd;
disp(['Investim_essnto em Imoveis/PIBd = ' num2str(Ih_PIBd)])
disp(' ')
K_PIBd = x(20)/PIBd;
disp(['Capital/PIBd = ' num2str(K_PIBd)])
disp(' ')
H_PIBd = (x(10) + x(11) + x(12))/PIBd;
disp(['Imóveis/PIBd = ' num2str(H_PIBd)])
disp(' ')
Emp_PIBd = (x(13) + x(14) + x(15))/PIBd;
disp(['(b_ih + b_iw + b_e)/PIBd = ' num2str(Emp_PIBd)])
disp(' ')
Dep_PIBd = (x(13) + x(14) + x(15) - x(23))/PIBd;
disp(['d/PIBd = ' num2str(Dep_PIBd)])
disp(' ')
disp(['b_ih/(b_ih + b_iw + b_e)(%): ',num2str((x(13)/(x(13) + x(14) + x(15)))*100)]);
disp(' ')
disp(['b_iw/(b_ih + b_iw + b_e)(%): ',num2str((x(14)/(x(13) + x(14) + x(15)))*100)]);
disp(' ')
disp(['b_e/(b_ih + b_iw + b_e)(%): ',num2str((x(15)/(x(13) + x(14) + x(15)))*100)]);
disp(' ')
r_d_anual = ((1 + r_d)^4 - 1);
disp(['r_d (%annual): ',num2str(r_d_anual)]);
disp(' ')
r_pol_anual = ((1 + r_pol)^4 - 1);
disp(['r_pol (%annual): ',num2str(r_pol_anual)]);
disp(' ')
r_bih_anual = ((1 + r_bih)^4 - 1);
disp(['r_bih (%annual): ',num2str(r_bih_anual)]);
disp(' ')
r_biw_anual = ((1 + r_biw)^4 - 1);
disp(['r_biw (%annual): ',num2str(r_biw_anual)]);
disp(' ')
r_be_anual = ((1 + r_be)^4 - 1);
disp(['r_be (%annual): ',num2str(r_be_anual)]);

end