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

NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
    paramname = M_.param_names{ii};
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;



A = 1;
KQ = 1;
ups = 1;

q = 1;

R = (1/betta-1)/(1-tau_r)+1;

% 目标
lev = 1.3663; % 相当于57.74负债率
F_cdf = 0.03; % 违约率
RkR = 1 + 0.02/4;% 溢价率


options = optimoptions('fsolve', 'OptimalityTolerance', 1e-7); % set options for numerical solver
x0 = [0.95, 0.3, 0.49, 0.008, 0.30, 0.1, 0.57, 1];
fun = @(var) myfun_findss(var,F_cdf,xi,lev,R,RkR);
var = fsolve(fun,x0,options);
gam        = var(1); 
sig_ome    = var(2); 
Gamma_cdf  = var(3);
G_cdf      = var(4);
F_pdf      = var(5);
mu         = var(6);
omega      = exp(var(7));
Rk         = var(8);

sig_ome_ss = sig_ome;

rk = Rk - (1 - delta);

kn = (rk/alphha)^(1/(alphha-1));

w = (1-alphha)*(kn)^(alphha); 

yn = kn^alphha;
cen = (1-gam)*(1-xi)*(1-Gamma_cdf)*Rk*kn;
xn = delta*kn;
cn = (1-gy_ss)*yn - cen - xn - mu*G_cdf*Rk*kn;
tlambda = (1-betta*h)*((1-h)*cn)^(-chi_C)/(1+tau_c);
Nd = ((1-tau_n)*tlambda*w/psi_N)^(1/(chi_N+chi_C));

k = kn*Nd;
y = yn*Nd;
Ce = cen*Nd;
x = xn*Nd;
C = cn*Nd;
lambda = tlambda/Nd^chi_C;

Ne = (gam + xi)*(1-Gamma_cdf)*Rk*k;
Loan = lev*Ne;

G = gy_ss*y;
T = G;
G_ss = G;


YY = 1;
KK = 1;
XX = 1;
CC = 1;
GG = 1;



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