function [ys,check] = NKopen_steadystate(ys,exe)
global M_ lgy_

if isfield(M_,'param_nbr') == 1
NumberOfParameters = M_.param_nbr;
for i = 1:NumberOfParameters
  paramname = deblank(M_.param_names(i,:));
  eval([ paramname ' = M_.params(' int2str(i) ');']);
end
check = 0;
end

options=optimset();
PI=PIbar;
u=1;
q=1;
d=1;
phi=1;
m=0;
zprime=exp(LambdaYd);
muprime=exp(Lambdamu);
Aprime=exp(LambdaA);
Lambdax=zprime;
r=(1-betta/zprime/muprime*(1-delta))/(betta/zprime/muprime);
gammma1=r;
R=PI*zprime/betta;
PIstar=((1-thetap*PI^(-(1-epsilon)*(1-chi)))/(1-thetap))^(1/(1-epsilon));
mc=(epsilon-1)/epsilon*(1-betta*thetap*PI^((1-chi)*epsilon))/(1-betta*thetap*PI^(-(1-epsilon)*(1-chi)))*PIstar;
PIstarw=((1-thetaw*PI^(-(1-chiw)*(1-eta))*zprime^(-(1-eta)))/(1-thetaw))^(1/(1-eta));
w=(1-alppha)*(mc*(alppha/r)^alppha)^(1/(1-alppha));
wstar=w*PIstarw;
vp=(1-thetap)/(1-thetap*PI^((1-chi)*epsilon))*PIstar^(-epsilon);
vw=(1-thetaw)/(1-thetaw*PI^((1-chiw)*eta)*zprime^eta)*PIstarw^(-eta);
tempvaromega=alppha/(1-alppha)*w/r*zprime*muprime;
ld=fsolve(@(ld)(1-betta*thetaw*zprime^(eta*(1+gammma))*PI^(eta*(1-chiw)*(1+gammma)))/(1-betta*thetaw*zprime^(eta-1)*PI^(-(1-chiw)*(1-eta)))...
-(psi*PIstarw^(-eta*gammma)*ld^gammma)/((eta-1)/eta*wstar*(1-h*betta*zprime^(-1))*(1-h/zprime)^(-1)*((Aprime/zprime*vp^(-1)*tempvaromega^alppha-(zprime*muprime-(1-delta))/(zprime*muprime)*tempvaromega)*ld-vp^(-1)*Phi)^(-1)),0.25,options);
l=vw*ld;
k=tempvaromega*ld;
x=(zprime*muprime-(1-delta))/(zprime*muprime)*k;
yd=(Aprime/zprime*k^alppha*ld^(1-alppha)-Phi)/vp;
c=(Aprime/zprime*vp^(-1)*tempvaromega^alppha-(zprime*muprime-(1-delta))/(zprime*muprime)*tempvaromega)*ld-vp^(-1)*Phi;
lambda=(1-h*betta*zprime^(-1))*(1-h/zprime)^(-1)/c;
f=(1-betta*thetaw*zprime^(eta-1)*PI^(-(1-chiw)*(1-eta)))^(-1)*(eta-1)/eta*wstar*PIstarw^(-eta)*lambda*ld;
g1=(1-betta*thetap*PI^((1-chi)*epsilon))^(-1)*lambda*mc*yd;
g2=epsilon/(epsilon-1)*g1;
Rbar=R;

% d*(c-h*c*zprime^(-1))^(-1)-h*betta*d*(c*zprime-h*c)^(-1)=lambda;
% 
% ((1-h/zprime)^(-1)-h*betta*1/(zprime-h))/c
for iter = 1:length(M_.params)
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

if isfield(M_,'param_nbr') == 1

if isfield(M_,'orig_endo_nbr') == 1
NumberOfEndogenousVariables = M_.orig_endo_nbr;
else
NumberOfEndogenousVariables = M_.endo_nbr;
end
ys = zeros(NumberOfEndogenousVariables,1);
for i = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(i,:));
  eval(['ys(' int2str(i) ') = ' varname ';']);
end
else
ys=zeros(length(lgy_),1);
for i = 1:length(lgy_)
    ys(i) = eval(lgy_(i,:));
end
check = 0;
end