function [ys,params,check] = mod_forum_steadystate(ys,exo,M_,options_)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% 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

omega        =   omega_ss;                                    
r            =   1/betta;                  
pi           =   1;                         
mc           =   ((psi1-1)/psi1)*(1/(1-omega)); 
alppha_n     =   1; 
Q = betta;

v=0;


Delta =   16.000000000000004;
B = 2.700000000000009;
%%%%%%%%%%%%%%%%% Solving for the subsystem of equations %%%%%%%%%%%%%%%%
x0 = [0.200000000000000; 0.900000000000000; 1.125000000000000;  0.225000000000000;    0.150000000000000;    0.650000000000000; 1; 0.027692307692308];  % guesses for Srf1, Sro1, S, Np
options = optimoptions('fsolve','TolFun',1e-15, 'MaxIter', 1e8,'Display','off','MaxFunctionEvaluations',500);
%x = fsolve(@root4_copy, x0, options, mc, theta1, v, alppha_n, sigma, chi_1, chi_2, chi_3, alppha_r, Tro1, Trf1, A , Delta)
x = fsolve(@root4_copy, x0, options, mc, theta1, v, alppha_n, sigma, chi_1, chi_2, chi_3, alppha_r, Tro1, Trf1, A , r, B, pi)
                                 



%%%%%%%%%% Initialising all the variables obtained with fsolve %%%%%%%%%%%%%%%
% x(1) = Np; x(2) = Srf1; x(3)= S;  x(4) = Sro1;  x(5) = Nro1;  x(6) = Nrf1
Np = x(1);
Srf1 = x(2);
S = x(3);
Sro1 = x(4);
Nro1 = x(5);
Nrf1 = x(6);
alppha_n = x(7);
t = x(8);




wf = (((1+alppha_n)*chi_1*(Nrf1)*((S)^(sigma)))/(1-t));
wo = (((1+alppha_n)*chi_2*(Nro1)*((S)^(sigma)))/(1-t));
wI = (((1+alppha_n)*chi_3*(Np)*((S)^(sigma)))/(1-t));

t=(v+(r*B)-B*pi)/(wf*Nrf1+wo*Nro1+wI*Np);

pI = wI/A;

Sirf1 = (theta1*(1-alppha_r)*mc*Srf1)/pI;

Siro1 = (theta1*(1-alppha_r)*mc*Sro1)/((1-v)*pI);

Si = Sirf1 + Siro1;

%Tro1=((mc*theta1)/((1-v)*((wo/alppha_r)^alppha_r)*((pI/(1-alppha_r))^(1-alppha_r))*((Sro1^((1-theta1)/theta1)))))^(1/alppha_r);
%Above equation just for Ramsey, not needed

% B = 2.700000000000009;
% 
% t=(v+(r*B)-B*pi)/(wf*Nrf1+wo*Nro1+wI*Np);
%B = (t*(wf*Nrf1+wo*Nro1+wI*Np)-v)/((r/pi)-1);
B_ss = B;

iota = t/((B/B_ss)^(xi));

r_ss         =   r;                                                                           
S_ss         =   S;                                                                           
Trf1_ss      =   Trf1;      

%A = Si/Np;   % just to make Ramsey work (not needed)

%% end own model equations


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
