function [ys,params,check] = model_dynamic_steadystate(ys,exo,M_,options_)

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%% Enter model equations here

if length(exo) >= 1
    flow_L = exo(1);
    flow_H = exo(2);
end

F_L_t2 = F0_L_t2 ;
F_H_t2 = F0_H_t2 + flow_H;
F_L_t1 = F0_L_t1 + flow_L;
F_H_t1 = F0_H_t1;
F = F_L_t2 + F_H_t2 + F_L_t1 + F_H_t1;

(F_H_t2 + F_L_t2)/F

r=(1 / bbeta) - 1;
rr =(1 - bbeta);

load('stst_guess.mat')
run solve_SS_HL_t1t2.m

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

u_H_t2 = (F_H_t2 * d_H)/(d_H + f_H);
n_H_t2 =  F_H_t2 - u_H_t2;
u_L_t2 = (F_L_t2 * d_L)/(d_L + f_L);
n_L_t2 =  F_L_t2 - u_L_t2;

u_H_t1 = (F_H_t1 * d_H)/(d_H + hH_t1 * f_H);
n_H_t1 =  F_H_t1 - u_H_t1;
u_L_t1 = (F_L_t1 * d_L)/(d_L + hL_t1 * f_L);
n_L_t1 =  F_L_t1 - u_L_t1;

u_L = u_L_t2 + u_L_t1;
u_H = u_H_t2 + u_H_t1;
v_H = u_H * theta_H;
v_L = u_L * theta_L;

q_H = f_H / theta_H;
q_L = f_L / theta_L;

Aj_L  = rr + bbeta * d_L;
Aj_H  = rr + bbeta * d_H;

l_L = Z_L * (n_L_t2 + n_L_t1 * lambda_L);
l_H = Z_H * ((chi * ((n_H_t1)^((epsilon - 1) / epsilon)) + (1 - chi) * ((n_H_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)

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;
i=r;

T = ((n_H_t1 * w_H_t1) * t + ...
     (n_H_t2 * w_H_t2) * t + ...
     (n_L_t1 * w_L_t1) * t + ...
     (n_L_t2 * w_L_t2) * t - ...
     (u_H_t1 * z_H_t1 + u_H_t2 * z_H_t2 + ...
      u_L_t1 * z_L_t1 + u_L_t2 * z_L_t2)) / F;
B_H_t1 = 0;
B_L_t1 = 0;
B_H_t2 = 0;
B_L_t2 = 0;
B = B_L_t1 + B_H_t1 + B_L_t2 + B_H_t2;

C_H_t2 = w_H_t2 * n_H_t2 * (1 - t) + (z_H_t2 * u_H_t2) + (F_H_t2 * T);
C_L_t1 = w_L_t1 * n_L_t1 * (1 - t) + (z_L_t1 * u_L_t1) + (F_L_t1 * T);
C_L_t2 = w_L_t2 * n_L_t2 * (1 - t) + (z_L_t2 * u_L_t2) + (F_L_t2 * T);

I = ddd * K;
C = Y - (I + (v_L * cL + v_H * cH));
C_H_t1 = C - C_H_t2 - C_L_t1 - C_L_t2;

syms W_H_t2_DE W_H_t1_DE W_L_t2_DE W_L_t1_DE U_H_t2_DE U_H_t1_DE U_L_t2_DE U_L_t1_DE

eqW_SS_1 = W_H_t1_DE == ((w_H_t1 * (1 - t)) * (1 / C_H_t1)) + bbeta * ((1 - d_H) * W_H_t1_DE + d_H * U_H_t1_DE);
eqW_SS_2 = W_H_t2_DE == ((w_H_t2 * (1 - t)) * (1 / C_H_t2)) + bbeta * ((1 - d_H) * W_H_t2_DE + d_H * U_H_t2_DE);
eqW_SS_3 = W_L_t1_DE == ((w_L_t1 * (1 - t)) * (1 / C_L_t1)) + bbeta * ((1 - d_L) * W_L_t1_DE + d_L * U_L_t1_DE);
eqW_SS_4 = W_L_t2_DE == ((w_L_t2 * (1 - t)) * (1 / C_L_t2)) + bbeta * ((1 - d_L) * W_L_t2_DE + d_L * U_L_t2_DE);
eqU_SS_1 = U_H_t1_DE == (z_H_t1 / C_H_t1) + bbeta * (f_H * hH_t1 * W_H_t1_DE + (1 - hH_t1 * f_H) * U_H_t1_DE);
eqU_SS_2 = U_H_t2_DE == (z_H_t2 / C_H_t2) + bbeta * (f_H * W_H_t2_DE + (1 - f_H) * U_H_t2_DE);
eqU_SS_3 = U_L_t1_DE == (z_L_t1 / C_L_t1) + bbeta * (f_L * hL_t1 * W_L_t1_DE + (1 - hL_t1 * f_L) * U_L_t1_DE);
eqU_SS_4 = U_L_t2_DE == (z_L_t2 / C_L_t2) + bbeta * (f_L * W_L_t2_DE + (1 - f_L) * U_L_t2_DE);

eqsUW_SS=[eqW_SS_1, eqW_SS_2, eqW_SS_3, eqW_SS_4, ...
          eqU_SS_1, eqU_SS_2, eqU_SS_3, eqU_SS_4];

varsUW_SS=[W_L_t2_DE, W_H_t2_DE, W_L_t1_DE, W_H_t1_DE, U_L_t1_DE, U_L_t2_DE, U_H_t1_DE, U_H_t2_DE];
          [W_L_t2, W_H_t2, W_L_t1, W_H_t1, U_L_t1, U_L_t2, U_H_t1, U_H_t2]=vpasolve(eqsUW_SS, varsUW_SS);

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
 