function [ys,check] = my_steadystate(junk,ys)
BETA      =   get_param_by_name('BETA');
BETA1     =   get_param_by_name('BETA1');
M         =   get_param_by_name('M');
MUC       =   get_param_by_name('MUC');
MUH       =   get_param_by_name('MUH');
DKC       =   get_param_by_name('DKC');
DKH       =   get_param_by_name('DKH');
DH        =   get_param_by_name('DH');
ETA       =   get_param_by_name('ETA');
ETA1      =   get_param_by_name('ETA1');
ALPHA     =   get_param_by_name('ALPHA');
NU        =   get_param_by_name('NU');
NU1       =   get_param_by_name('NU1');
KAPPA     =   get_param_by_name('KAPPA');
MUBB      =   get_param_by_name('MUBB');
sigma   =    get_param_by_name('sigma');
omega   =    get_param_by_name('omega');
rho   =    get_param_by_name('rho');
xi_p  =    get_param_by_name('xi_p');

% equations DEFINING CAPITAL TO PRODUCTION RATIOS (1) to (6) in "Steady state"
check=0;

%%

% Interest rates
r = 1 / BETA ;
rkc = r - (1-DKC) ;
rkh = r - (1-DKH) ;
r_f =  r ;
r_F = r-1;
r1 = r-1;

ZETA0 = BETA*MUC/(1-BETA*(1-DKC)) ;
ZETA1 = BETA*MUH/(1-BETA*(1-DKH));
ZETA2 = 1/(1-BETA*(1-DH)) ;
ZETA3 = 1/(1-BETA1*(1-DH)-(BETA-BETA1)*M) ;
ZETA4 = (r-1)*M/r ;

DH1 = 1 - (1-DH) ;
DKC1 = 1 - (1-DKC) ; 
DKH1 = 1 - (1-DKH) ;

CHI1 = 1+DH1*ZETA2*(1-r1*ZETA1-KAPPA-ALPHA*(1-MUH-KAPPA-MUBB)) ;
CHI2 = (r1*ZETA1+KAPPA+ALPHA*(1-MUH-KAPPA-MUBB))*DH1*ZETA3+ZETA4*ZETA3 ;
CHI3 = (r1*ZETA0+ALPHA*(1-MUC));
CHI4 = 1+DH1*ZETA3*(1-(1-ALPHA)*(1-MUH-KAPPA-MUBB))+ZETA4*ZETA3 ;
CHI5 = (1-ALPHA)*(1-MUH-KAPPA-MUBB)*DH1*ZETA2 ;
CHI6 = (1-ALPHA)*(1-MUC) ;


CY = (CHI3*CHI4+CHI2*CHI6)/(CHI1*CHI4-CHI2*CHI5) ;
CYPRIME = (CHI1*CHI6+CHI3*CHI5)/(CHI1*CHI4-CHI2*CHI5) ;
QIY = DH1*ZETA2*CY + DH1*ZETA3*CYPRIME ;

RATION = (1-MUH-KAPPA-MUBB)/(1-MUC)*QIY  ;
NHNC = RATION^(1/(1-NU)) ;
NHNC1 = RATION^(1/(1-NU1)) ;

nc = ( ((1-MUC)*ALPHA/CY)/(1+RATION)^((ETA+NU)/(1-NU)) )^(1/(1+ETA)) ;
nh = NHNC*nc ;

nc1 = ( ((1-MUC)*(1-ALPHA)/CYPRIME)/(1+RATION)^((ETA1+NU1)/(1-NU1)) )^(1/(1+ETA1)) ;
nh1 = NHNC1*nc1 ;

Y = (nc^ALPHA)*(nc1^(1-ALPHA)) *  ZETA0^(MUC/(1-MUC))  ;
y = Y;
y_d = omega*y;
y_m = (1-omega)*y;
y_x = (1-omega)*y;

QI = QIY*Y ;
kh = ZETA1*QI ;
BRI = MUBB*QI ;
I = (nh^ALPHA*nh1^(1-ALPHA))^(1-MUH-MUBB-KAPPA)*kh^MUH*BRI^MUBB;
q = QIY*Y / I ;

kc = ZETA0*Y ;

c = CY*Y ;
c1 = CYPRIME*Y ;

h = ZETA2*c/q ;
h1 = ZETA3*c1/q ;

b = M*q*h1/r ;

lm = (1-BETA1/BETA)/c1 ;

uc = 1/c;
uc1 = 1/c1;

K_d = uc*(1-1/sigma)*y/(1-BETA*xi_p);
F_d = uc*y/(1-BETA*xi_p);

d_p = 1 ;
bf  = 0;
mc  =  1;
pi = 1;
pi_s = 1;
pi_d = 1;
pi_m = 1;
pi_star = 1;
tau = 0;
cm = 1;
T = 1;

wc = ALPHA*(1-MUC)*Y/nc ;
wh = ALPHA*(1-MUH-KAPPA-MUBB)*QI/nh ;
wc1 = (1-ALPHA)*(1-MUC)*Y/nc1 ;
wh1 = (1-ALPHA)*(1-MUH-KAPPA-MUBB)*QI/nh1 ;

xxx=[
    c
    c1
    1
    b
    mc
    d_p
    h
    h1
    I
    kc
    kh
    lm
    nc
    nc1
    nh
    nh1
    q
    r
    rkc
    rkh
    r_f
    exp(r_F)
    uc
    uc1
    wc
    wc1
    wh
    wh1
    Y
    y
    y_d
    y_x
    y_m
    pi
    pi_s
    pi_d
    pi_m
    pi_star
    K_d
    F_d
    1
    exp(1)
    exp(1)];

if (any(xxx<=0))
  % if the steady state does not exist, make penalty from negative
  % numbers and return check=-1, this is handled in resol.m, and the
  % penalty is put to info(2)
  penalty = sum(xxx(find(xxx <= 0)).^2);
  ys = zeros(size(xxx));
  ys(1) = penalty;
  check = -1;
else
  logxxx = log(xxx) ;
  ys = logxxx ;
end
