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


 global M_
 
 BETAH = get_param_by_name('BETAH') ;
 BETAB = get_param_by_name('BETAB') ;
 BETAE = get_param_by_name('BETAE') ;
 BETAS = get_param_by_name('BETAS') ;
 ALPHA = get_param_by_name('ALPHA') ;
 MS = get_param_by_name('MS') ;
 MH = get_param_by_name('MH') ;
 MK = get_param_by_name('MK') ;
 MN = get_param_by_name('MN') ;
 GAMMAE = get_param_by_name('GAMMAE') ;
 GAMMAS = get_param_by_name('GAMMAS') ;
 DKH = get_param_by_name('DKH') ;
NU =   get_param_by_name('NU');
RHOD =   get_param_by_name('RHOD');
RHOE =   get_param_by_name('RHOE');
RHOS =   get_param_by_name('RHOS');
SIGMA =   get_param_by_name('SIGMA');
SIGMAB =   get_param_by_name('SIGMAB');
SIGMAE =   get_param_by_name('SIGMAE');
SIGMAH =   get_param_by_name('SIGMAH');
GAMMA = get_param_by_name('GAMMA');
DELTA = get_param_by_name('DELTA') ;
OMEGA = get_param_by_name('OMEGA');
PSI = get_param_by_name('PSI');
THETA = get_param_by_name('THETA');

check=0;
BETA_H=BETAH * (1-THETA*(1-exp(1-GAMMA)*(1-DELTA)))^((1-1/(1/(1-(1+OMEGA)*(1-PSI))))/((1-GAMMA)*(1+OMEGA)));
BETA_S=BETAS * (1-THETA*(1-exp(1-GAMMA)*(1-DELTA)))^((1-1/(1/(1-(1+OMEGA)*(1-PSI))))/((1-GAMMA)*(1+OMEGA)));
rh = 1/BETA_H; 
rkh = rh - (1-DKH);

lambdab = (1-BETAB*rh)/(1-BETAB*RHOD) ;

re = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAE)*lambdab) ; 
rs = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAS)*lambdab) ; 

lambdae = (1-BETAE*re)/(1-BETAE*RHOE) ;
lambdas = (1-BETA_S*rs)/(1-BETA_S*RHOS) ;
rv = 1/BETAE*(1-lambdae*MH/re*(1-RHOE)) - 1   ;

oo1 = 1/(1-BETA_H) ;
oo2 = 1/(1-BETA_S-lambdas*(1-RHOS)*MS/rs);
oo3 = 1/(1+(1-1/rs)*MS*oo2);
oo4 = GAMMAE*(rh-1)*(MH*NU/(re*rv)-MN*(1-ALPHA-NU)/(1+MN*lambdae)) ;
oo5 = GAMMAS*MS/rs*(oo2*(rh-1))*oo3 ;
oo6 = rkh - DKH ;
oo7 = (1-ALPHA-NU)*(1-SIGMA)/((1+MN*lambdae)) ;
oo8 = ALPHA/(rkh) ;


ns = 1/(1+oo3*OMEGA) ;

zeta1 = oo6*oo8/oo7 + 1 + oo4/oo7 + oo5*SIGMA/(1-SIGMA)  ;

nh = 1/(1+OMEGA*zeta1) ;

ctoys = oo3*(1-ALPHA-NU)/(1+MN*lambdae)*SIGMA ;
ctoyh = zeta1*oo7 ;

iktoyh = DKH*ALPHA/(rkh) ;

denh = oo1*ctoyh + oo2*ctoys + NU/(rv) ;
hh = oo1*ctoyh / denh ;
hs = oo2*ctoys / denh ;
he = (NU/(rv)) / denh ;

y = (ALPHA/(rkh))^(ALPHA/(1-ALPHA))*...
    he^(NU/(1-ALPHA))*...
    (nh^(1-SIGMA)*ns^SIGMA)^((1-ALPHA-NU)/(1-ALPHA)) ;

kh = ALPHA*y/(rkh) ;


q = NU*y/(rv)/he ;

loe = MH*q/re*he - MN*(1-ALPHA-NU)/(1+MN*lambdae)*y ;

wh = ((1-SIGMA)*(1-ALPHA-NU)/(1+MN*lambdae))*y/nh ;
ws = (SIGMA*(1-ALPHA-NU)/(1+MN*lambdae))*y/ns;

los = MS/rs*oo2*oo3*ws*ns ;

d = GAMMAE*loe + GAMMAS*los ;



ch = (rkh-DKH)*kh + (rh-1)*d + wh*nh   ;

cb = (re-1)*loe + (rs-1)*los - (rh-1)*d ;

ce = y -  (re-1)*loe - wh*nh - ws*ns - rkh*kh ;
cs = oo3*ws*ns ;


dp = 1;

ctot = ch + cs + ce + cb;
itot = DELTA*(kh) ;

kb = loe + los - d;

spread_re_rh = 4*(re-rh);
spread_rs_rh = 4*(rs-rh);

ktot = kh;
ntot =  nh^(1-SIGMA)*ns^SIGMA  ;


ytot = itot+ctot;


acdb=0;
acdh=0;
ache=0;
achh=0;
achs=0;
acke=0;
ackh=0;
acloeb=0;
aclosb=0;
acloee=0;
acloss=0;

Y_SS = (y);
ITOT_SS = (itot);
CTOT_SS = (ctot);
Q_SS = (q);
LTOT_SS = loe+los;
NTOT_SS = ntot;
KH_SS = kh;
D_SS = d;
LOE_SS = loe;
LOS_SS = los;

intshare = 1;

loetoy=loe/4/y;
lostoy=los/4/y;
lotoy=(loe+los)/4/y;

cbtoy = cb/y ;
cetoy = ce/y ;
chtoy = ch/y ;
cstoy = cs/y ;


kab=1-(d/(loe+los));

data_CC=0;
data_KAB=kab;
data_II=0;
data_LL=0;
data_LO=0;
data_LOE=0;
data_LOH=0;
data_QQ=0;
data_YY=0;
data_YL=0;
data_TFP=0;
data_NN=0;

zke=1;
zkh=1;
dkh=DKH;

CE_SS = ce;
RE_SS = re;
LE_SS = lambdae;
LB_SS = lambdab;
CB_SS = cb;

ucb = cb^-SIGMAB;
uce = ce^-SIGMAE;

a_be=0;
a_bh=0;
a_j=0;
a_k=0;
a_me=0;
a_mh=0;
a_p=0;
a_z=0;
theta =0;

e_a_be=0;
e_a_bh=0;


xxx = [ ...
a_be
a_bh
a_j
a_k
a_me
a_mh
a_p
a_z
acdb
acdh
ackh
acloeb
acloee
aclosb
acloss
cb
ce
ch
cs
ctot
d
data_CC
data_II
data_LO
data_LOE
data_LOH
data_NN
data_QQ
data_TFP
data_YY
e_a_be
e_a_bh
he
hh
hs
itot
kh
ktot
lambdab
lambdae
lambdas
loe
los
nh
ns
ntot
q
re
rh
rkh
rs
rv
spread_re_rh
spread_rs_rh
ucb
uce
wh
ws
y 
BETA_H
BETA_S
theta] ;



ys = (xxx) ;

