function [ys,params,check] = lianxi4_steadystate(ys,exo,M_,options_)
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] 0 if steady state computation worked and to
%                        1 of not (allows to impose restrictions on parameters)

%% Step 0: initialize indicator and set options for numerical solver
check = 0;
options = optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

%% Step: read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
    paramname = deblank(M_.param_names{ii,:});
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
%% 

  mu = muss;
  Rak = 1/betta - 1 + del_ak;
  Raz = 1/betta - 1 + del_az;
  Yag = 1;
  Yang = 1;
  Na = 1/3;
  Ya = ((xaj)^(1/omegaa)+(1-xaj)^(1/omegaa))^(omegaa/(omegaa-1));
  Pag = (xaj * Ya)^(1/omegaa);
  Pang = ((1-xaj) * Ya)^(1/omegaa);
  Kag = ((omegaa-1)/omegaa) * Pag * alpag / Rak ;
  Kang = ((omegaa-1)/omegaa) * Pang * alpang / Rak ;
  %Vag Vang Nag Nang Zag Zang wa
  fun_a = @(x)Fun_a(x,e1,sig1,omegaa,Pag,alpag,e2,sig2,Pang,alpang,Raz,Na);
  x0=[0.3,0.3,0.1,0.2,0.4,0.5,1.5];
  options = optimoptions('fsolve', 'Display', 'none','TolX', 1e-8,'TolFun', 1e-8);
  x = fsolve(fun_a, x0, options);
  Vag = x(1);
  Vang = x(2);
  Nag = x(3);
  Nang = x(4);
  Zag = x(5);
  Zang = x(6);
  wa = x(7);
  Aag = Yag/(Kag^alpag * Vag^(1-alpag));
  Aang = Yang/(Kang^alpang * Vang^(1-alpang));
  Ka = Kag + Kang;
  Za = Zag + Zang;
  Iak = del_ak * Ka;
  Iaz = del_az * Za;
  G = mu * Yag;
  Ca = Ya - del_ak *  Ka -del_az *  Za - G;
  omega = wa/(Ca^(sig_ac) * Na^sig_an);

%% Step: Update parameters and variables
params=NaN(NumberOfParameters,1);
for iter = 1:M_.param_nbr %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


end
function F=Fun_a(x,e1,sig1,omegaa,Pag,alpag,e2,sig2,Pang,alpang,Raz,Na)
  Vag = x(1);
  Vang = x(2);
  Nag = x(3);
  Nang = x(4);
  Zag = x(5);
  Zang = x(6);
  wa = x(7);

F(1)=Vag - ((1 - e1) * Nag^((sig1 - 1)/sig1) + e1 * Zag^((sig1 - 1)/sig1))^(sig1/(sig1 - 1));
F(2)=Vang - ((1 - e2) * Nang^((sig2 - 1)/sig2) + e2 * Zang^((sig2 - 1)/sig2))^(sig2/(sig2 - 1));
F(3)= wa - ((omegaa-1)/omegaa) * Pag * (1-alpag)* (1-e1)  / Vag * (Vag/ Nag)^(1/sig1);
F(4) = Raz - ((omegaa-1)/omegaa) * Pag * (1-alpag) * e1 / Vag * (Vag/ Zag)^(1/sig1); 
F(5) = wa - ((omegaa-1)/omegaa) * Pang * (1-alpang) * (1-e2)  / Vang * (Vang/ Nang)^(1/sig2); 
F(6) = Raz - ((omegaa-1)/omegaa) * Pang * (1-alpang)* e2 / Vang * (Vang/ Zang)^(1/sig2); 
F(7) =  Nag + Nang - Na;
end