function [ys,check] = rbc_3_steadystate(ys,exe)
  global M_
  
  beta = M_.params(1);
  B=M_.params(2);
  delta = M_.params(3);
  theta=M_.params(4);
  alfa = M_.params(5);
  RHO_a = M_.params(6); 
  RHO_q = M_.params(7);
  RHO_c = M_.params(8);
  MU_a = M_.params(9);
  MU_q = M_.params(10);
  MU_c = M_.params(11);
  A_n = M_.params(12); 
   

  check = 0;
  

  da=exp(MU_a);
  MU_aa=(exp(MU_a)^(-1/(1+alfa)));
  MU_aaa=(exp(MU_a)^(-1/(1-alfa)));
  MU_aaaa=(((1-alfa)/A_n)^((1-alfa)/(1+alfa)))*(MU_c^((-1)/(1+alfa)))*((1-B/exp(MU_a))^((-1)/(1+alfa)));
  
  xx=(beta*alfa*MU_aaa*(exp(MU_a)-(1-delta))+(1-beta*MU_aaa*(1-delta))*(alfa+MU_q)*theta)/((1-beta*MU_aaa*(1-delta))*(alfa+MU_q)*theta-beta*alfa*MU_aaa*(exp(MU_a)-(1-delta))*theta);
  
  x=xx^(1/(1+theta));
  kk=MU_aaaa*(alfa+MU_q)*theta*(xx-1)/(1+theta)/xx/(exp(MU_a)-(1-delta));
  k=kk^(1+alfa);
  y=MU_aaaa*(kk^alfa);
  i=y*(alfa+MU_q)*(1-x^(-theta));
  c=MU_c;
  q=MU_q;
  x_a=0;
  di=da;
  dy=da;
  
 
  
  ys =[
da     
x     
k     
y     
i     
c     
q     
di        
x_a
dy
];