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




rhoa = M_.params(1);
rhox = M_.params(2);
rhoz = M_.params(3);
rhogg = M_.params(4);
thetac = M_.params(5);
thetal = M_.params(6);
beta = M_.params(7);
alfa = M_.params(8);
constelab = M_.params(9);
delta0 = M_.params(10);
delta1 = M_.params(11);
delta2= M_.params(12);
sigma= M_.params(13);
rhoxg= M_.params(14);
miua= M_.params(15);
GY= M_.params(16);
kai= M_.params(17);
kappa= M_.params(18);
miux = M_.params(19);
miuk = M_.params(20);
miuy= M_.params(21);
miug= M_.params(22);


check = 0;

v_u = 1;
v_h = ( (1-GY-(miuk-1+delta0)*alfa/delta1)/(((1-alfa)*(1-thetal)*(1-thetac*beta*miuy^(-sigma)))/(kai*(1-thetac/miuy)*(1-thetal*beta*miuy^(1-sigma))))+ 1)^(-1);
v_l = 1-v_h;
v_miux = miux;
v_miuk = miuk;
v_miuy = miuy;
v_miua = miua;
v_k = ((((1/(beta*miua*miuy^(-sigma)))+delta0-1)/alfa)^(1/(alfa-1)))*miuk*v_h; 
v_y = (v_k/miuk)^alfa*(v_h)^(1-alfa);
v_xg = (miuy)^(-1/(1-rhoxg*0.99));
v_g = GY*v_y/v_xg;
v_i = v_k - (1-delta0)*v_k/miuk;
v_z = 1;
v_c = v_y - v_i - v_g*v_xg;
v_lamda = (((1-thetac/v_miuy)*v_c)*((1-thetal)*v_l)^kai)^(-sigma)*((1-thetal)*v_l)^kai -thetac*beta*(((v_miuy-thetac)*v_c)*((1-thetal)*v_l)^kai)^(-sigma)*((1-thetal)*v_l)^kai;
v_q = 1;
v_an1 = 0;
v_an2 = 0;
v_an3 = 0;
v_an2lag1 = 0;
v_an3lag1 = 0;
v_an3lag2 = 0;
dy = v_miuy;
dc = v_miuy;
dinve = v_miuk;
eks = v_miuk;
eks1 = v_miuk;
es = 0;
es1 = 0;
esprime =0;
esprime1 =0;


ys = [
v_k 
v_u 
v_miuk 
v_i 
v_y 
v_z 
v_h 
v_lamda 
v_l 
v_c 
v_miuy 
v_g 
v_q 
v_miua 
v_miux 
v_xg      
v_an1 
v_an2 
v_an3 
v_an2lag1 
v_an3lag1 
v_an3lag2 
dy 
dc 
dinve	  	   
eks 
eks1 
es 
es1 
esprime 
esprime1 ];