function [ys,params,check] = Policy_Rep_steadystate(ys,exo,M_,options_)
  
NumberOfParameters = M_.param_nbr;
[alpha,betta,gamma_d,delta,Delta,eta,theta,theta_g,theta_b,theta_B,theta_I,theta_P,theta1_g,theta2_g,theta1_b,theta2_b,kappa,lambda,lambda_g,lambda_b,pi_bar,rho_c,rho_epsiA,sigma,epsilonA_bar,phi_g,phi_b,phi_pi,phi_y,omega,a,b,Emi_Star,g_y,hab,i_bar,L,v1,v2]=deal(NaN);
for ii = 1:NumberOfParameters
paramname = M_.param_names{ii};
eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

check=0;

options=optimset('Display','off');

R=1/betta;
pi=1;
pi_g=1;
pi_b=1;
Q_g=1;
Q_b=1;
P=1;
i=i_bar;
epsilonA=epsilonA_bar;
Omega=(1-theta_B)/(1-Delta*betta*(1-omega));
Gamma_B=(Omega+theta_B-1)/theta_B;
nu=1-Delta*Omega/Gamma_B;
RP_g=nu*lambda*lambda_g/(Delta*betta*Omega);
RP_b=nu*lambda*lambda_b/(Delta*betta*Omega);
R_g=RP_g+R;
R_b=RP_b+R;
RK_g=R_g-1+delta;
RK_b=R_b-1+delta;

x0=[0.12;0.23;0.4;0.4;1.5;1.5;22;22;1;1;3;3;1;1;0.9;0.1;2;1700;0.2;0.1;0.5];
[y,fval,exitflag]=fsolve(@(x) Policy_Rep_steady_state_helper(x,phi_g,phi_b,a,b,alpha,theta1_g,theta2_g,theta1_b,theta2_b,v1,v2,Emi_Star,gamma_d,kappa,L,RK_g,RK_b,betta),x0,options);
if exitflag <1
    check=1;
    return;
end
Emi_g=y(1);
Emi_b=y(2);
mu_g=y(3);
mu_b=y(4);
Y_g=y(5);
Y_b=y(6);
K_g=y(7);
K_b=y(8);
W_g=y(11);
W_b=y(12);
MC_g=y(13);
MC_b=y(14);
da=y(15);
tau=y(16);
T=y(17);
X=y(18);
Emi=y(19);
SCC=y(20);
ksi=y(21);

P_g=theta_g*MC_g/(theta_g-1);
P_b=theta_b*MC_b/(theta_b-1);
I_g=delta*K_g;
I_b=delta*K_b;
S_g=K_g;
S_b=K_b;
N=theta_B*(RP_g*S_g+RP_b*S_b)/(1-theta_B*R-omega);
V=Gamma_B*N;
I=kappa*I_g+(1-kappa)*I_b;
K=kappa*K_g+(1-kappa)*K_b;
Z_g=theta1_g*mu_g^theta2_g*Y_g;
Z_b=theta1_b*mu_b^theta2_b*Y_b;
Z=kappa*Z_g+(1-kappa)*Z_b;
Y=(kappa*Y_g^(1-1/theta)+(1-kappa)*Y_b^(1-1/theta))^(1/(1-1/theta));
G=g_y*Y;
Tax=G-tau*Emi;
C=Y-I-G-Z;
rho=(1-betta*hab)*(C-hab*C)^-sigma;

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params)
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr;
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end