function [ys,check] = model_steadystate(ys,exo)
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





%%%% Model equations to be entered here




lambdass=2;

lambda=lambdass;

z_bar=lambda^(1/eta);

Z=(eta/(eta-1)*z_bar)^alppha;

S=0;



CW=(alppha/(1/betta-(1-delta)))^(alppha/(1-alppha))*Z^(1/(1-alppha))*(1-alppha)*(h1^(1/alppha)*(a1-a0)^((alppha-1)/alppha)+h2^(1/alppha)*(a2-a1)^((alppha-1)/alppha))^(alppha/(1-alppha));



LF_K =((1-alppha)*h1*Z/(a1-a0)*CW)^(1/alppha); %LF/K

LP_K = ((1-alppha)*h2*Z/(a2-a1)*CW)^(1/alppha); %LP/K

w1=(1-alppha)*Z*(LF_K)^(-alppha);

w2=(1-alppha)*Z*(LP_K)^(-alppha);



Y=Z*K^alppha*(h1*LF^(1-alppha)+h2*(LP)^(1-alppha)); % can not get analytical solution

K=alppha*Y/(1/betta-(1-delta));                     % can not get anlaytical solution



%LF = (((1-alppha)*h1*Z)/((a1-a0)*CW))^(1/alppha)*K;
%LP = (((1-alppha)*h2*Z)/((a2-a1)*CW))^(1/alppha)*K;


%K = betta * (alppha * Y + (1-delta)*K ); 

%w1=(1-alppha)*Z*K^alppha*LF^(-alppha);

%w2=(1-alppha)*Z*K^alppha*LP^(-alppha);

%CW=w1*h1*LF+w2*h2*LP;

CE=(1-betta)*(alppha*Y+(1-delta)*K);

C=CW+CE;

pi=alppha*(h1)^(1/alppha)*((1-alppha)/(w1*h1))^((1-alppha)/alppha)+alppha*h2*((1-alppha)/w2)^((1-alppha)/alppha);

r=pi*z_bar-delta;

I=delta*K;







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
end