function [ys,check]=replic_steadystate(junk,ys)


 global M_
 
 beta = get_param_by_name('beta') ;
 beta2 = get_param_by_name('beta2') ; 
 gamma = get_param_by_name('gamma') ;
 j = get_param_by_name('j') ;
 eta = get_param_by_name('eta') ;
 my = get_param_by_name('my') ;
 ypsilon = get_param_by_name('ypsilon') ;
 psi = get_param_by_name('psi') ;
 delta = get_param_by_name('delta') ;
 X = get_param_by_name('X') ;
 theta = get_param_by_name('theta') ;
 alfa = get_param_by_name('alfa') ;
 m = get_param_by_name('m') ;
 mss = get_param_by_name('mss') ;
 rhou = get_param_by_name('rhoj') ;
 rhoA = get_param_by_name('rhoA') ;
 sigmau = get_param_by_name('sigmau') ;
 sigmaj = get_param_by_name('sigmaj') ;
 sigmaA = get_param_by_name('sigmaA') ;
 sigmaR = get_param_by_name('sigmaR') ;
 

 check=0;
 
 % Calculating steady state of the model
z1 = (gamma*my) / ( X*(1-gamma*(1-delta)) ) ;
z2 = (gamma*ypsilon) /( X*(1-gamma-m*(beta-gamma)) );
s2=(1-alfa)*(1-my-ypsilon)/X; 
s1=(alfa*(1-my-ypsilon)+X-1)/X; % income shares
z3 =  j/(1-beta)  ;
z4 =  j/(1-beta2-mss*(beta-beta2)) ;
z5 =  (my+ypsilon)/X - delta*z1 - (1-beta)*m*z2  ;  
z6 = s2 / ( 1+(1-beta)*mss*z4 ) ;
z7 = s1 + (1-beta)*( m*z2 + mss*z4*z6 )   ; 

denh = z3*z7+z4*z6+z2 ;

Y = 1;

h = z2/denh ; 
h1 = z3*z7/denh ; 
h2 = z4*z6/denh ;
q = z2/h ; 

K = z1 ;  
c = z5 ;  
c2 = z6 ;  
c1 = z7  ; 

b = beta*m*q*h ;

r = 1/beta;
pi = 1;

b2 = beta*mss*q*h2 ;
lambda = (beta-gamma)/c ; lambda2 = (beta-beta2)/c2 ;
I = delta*K ;

b1 = -b2 - b;

L1 = (s1/c1)^(1/eta);
L2 = (s2/c2)^(1/eta);

w1 = s1/L1;
w2 = s2/L2;


a_j = 0;
a_u = 0;
a_a = 0;

 
xxx = [ ...
Y,
h,
h1,
h2,
q,
K,
c,
c1,
c2,
r,
pi,
b,
b1,
b2,
lambda,
lambda2,
I,
w1,
L1,
w2,
L2,
a_j,
a_a,
a_u ] ;



ys = (xxx) ;