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


 global M_
 
 BETAH = get_param_by_name('BETAH') ;
 BETAB = get_param_by_name('BETAB') ;
 BETAI = get_param_by_name('BETAI') ;
 BETAS = get_param_by_name('BETAS') ;
 BETAN = get_param_by_name('BETAN') ;
 ALPHA = get_param_by_name('ALPHA') ;
 MS = get_param_by_name('MS') ;
 MHS = get_param_by_name('MHS') ;
 MHN = get_param_by_name('MHN') ;
 MKS = get_param_by_name('MKS') ;
 MKN = get_param_by_name('MKN') ;
 MNS = get_param_by_name('MNS') ;
 MNN = get_param_by_name('MNN') ;
 GAMMAI = get_param_by_name('GAMMAI') ;
 GAMMAN = get_param_by_name('GAMMAN') ;
 GAMMAS = get_param_by_name('GAMMAS') ;
 JEI = get_param_by_name('JEI') ;
 DELTA = get_param_by_name('DELTA') ;
 TAUH = get_param_by_name('TAUH') ;
 TAUI = get_param_by_name('TAUI') ;

MIU    = get_param_by_name('MIU') ;
NU =   get_param_by_name('NU');
RHOD =   get_param_by_name('RHOD');
RHOI =   get_param_by_name('RHOI');
RHON =   get_param_by_name('RHON');
RHOS =   get_param_by_name('RHOS');
SIGMA =   get_param_by_name('SIGMA');

SIGMAB =   get_param_by_name('SIGMAB');
SIGMAI =   get_param_by_name('SIGMAI');
SIGMAH =   get_param_by_name('SIGMAH');
SIGMAN =   get_param_by_name('SIGMAN');
SIGMAS =   get_param_by_name('SIGMAS');
DELTAKH = get_param_by_name('DELTAKH');
DELTAKS = get_param_by_name('DELTAKS');
DELTAKN = get_param_by_name('DELTAKN');
MI = get_param_by_name('MI');
kappa = get_param_by_name('kappa');
Epsilon = get_param_by_name('Epsilon');
check=0;

rh = 1/BETAH; 
rkh = rh - (1-DELTAKH);

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

ri = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAI)*lambdab) ; 
rs = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAS)*lambdab) ; 
rn = (1/BETAB)*(1-((1-BETAB)*RHOD+(1-RHOD)*GAMMAN)*lambdab) ; 

lambdai = (1-BETAI*ri)/(1-BETAI*RHOI) ;
lambdan = (1-BETAN*rn)/(1-BETAN*RHON) ;
lambdas = (1-BETAS*rs)/(1-BETAS*RHOS) ;

rks = 1/BETAS*(1-lambdas*(1-RHOS)*MKS) - (1-DELTAKS) ;
rvs = 1/BETAS*(1-lambdas*MHS/rs*(1-RHOS)) - 1   ;

rkn = 1/BETAN*(1-lambdan*(1-RHON)*MKN) - (1-DELTAKN) ;
rvn = 1/BETAN*(1-lambdan*MHN/rn*(1-RHON)) - 1   ;

oo1 = JEI/(1-BETAH) ;
oo2 = JEI/(1-BETAI-lambdai*(1-RHOI)*MI/ri);
oo3 = 1/(1+(1-1/ri)*MI*oo2);
oo4 = GAMMAS*(rh-1)*(MHS*NU/(rs*rvs)+MKS*ALPHA/(rks*rs)-MNS*(1-ALPHA-NU)/(1+MNS*lambdas)) ;
oo41 = GAMMAN*(rh-1)*(MHN*NU/(rn*rvn)+MKN*ALPHA*(1-MIU)/(rkn*rn)-MNN*(1-ALPHA-NU)/(1+MNN*lambdan)) ;
oo5 = GAMMAI*MI/ri*(oo2*(rh-1))*oo3 ;
oo6 = rkh - DELTAKH ;
oo7 = (1-ALPHA-NU)*(1-SIGMA)/(kappa*(1+MNS*lambdas)) ;
oo71 = (1-ALPHA-NU)*(1-SIGMA)/((1-kappa)*(1+MNN*lambdan)) ;
oo8 = MIU*ALPHA/(rkh) ;


ni = 1/(1+oo3*TAUI) ;

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

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

denh = oo1*oo7*zeta1 + oo2*oo3*oo7*SIGMA/(1-SIGMA) + NU/(rvs) + NU*oo7/(rvn*oo71) ;

hh = oo1*zeta1*oo7 / denh ;
hi = oo2*oo3*oo7*SIGMA/(1-SIGMA) / denh ;
hs = (NU/(rvs)) / denh ;
hn = ((NU*oo7)/(rvn*oo71)) / denh ;

yn = (ALPHA*(1-MIU)/(rkn))^(ALPHA*(1-MIU)/(1-ALPHA))* (ALPHA*MIU/(rkh))^(ALPHA*MIU/(1-ALPHA))* hn^(NU/(1-ALPHA))* (((1-kappa)*nh)^(1-SIGMAH)*((1-kappa)*ni)^SIGMAI)^((1-ALPHA-NU)/(1-ALPHA)) ;

ys = Epsilon^(1/(1-ALPHA)) * (ALPHA/(rks))^(ALPHA/(1-ALPHA))* hs^(NU/(1-ALPHA))* ((kappa*nh)^(1-SIGMAH)*(kappa*ni)^SIGMAI)^((1-ALPHA-NU)/(1-ALPHA)) ;
 
kh = ALPHA*MIU*yn/(rkh) ;
ks = ALPHA*ys/(rks) ;
kn = ALPHA*(1-MIU)*yn/(rkn) ;
q = NU*ys/(rvs)/hs ;


los = MHS*q/rs*hs + MKS*ks - MNS*(1-ALPHA-NU)/(1+MNS*lambdas)*ys ;
lon = MHN*q/rn*hn + MKN*kn - MNN*(1-ALPHA-NU)/(1+MNN*lambdan)*yn ;
loi = MI*q*hi/ri ;

wh = ((1-SIGMA)*(1-ALPHA-NU)/(kappa*(1+MNS*lambdas)))*ys/nh ;
wi = (SIGMA*(1-ALPHA-NU)/(kappa*(1+MNS*lambdas)))*ys/ni;

d = GAMMAI*loi + GAMMAS*los + GAMMAN*lon ;
ch = (rkh-DELTAKH)*kh + (rh-1)*d + wh*nh ;
cb = ((ri-1)*loi + (rs-1)*los + (rn-1)*lon - (rh-1)*d);


mu = cb^-SIGMAB ;



cs = ys - DELTAKS*ks - (rs-1)*los - kappa*wh*nh - kappa*wi*ni ;
cn = yn - DELTAKN*kn - (rn-1)*lon - (1-kappa)*wh*nh - (1-kappa)*wi*ni -rkh*kh;
ci = oo3*wi*ni ;

dp = 1;

ctot = ch + cs + ci + cs + cb;
itot = DELTAKH*kh+DELTAKS*ks+DELTAKN*kn ;

kb = loi + los + lon - d;

spread_ri_rh = 4*(ri-rh);
spread_rn_rh = 4*(rn-rh);
spread_rs_rh = 4*(rs-rh);

ktot = ks+kh+kn;
ntot = nh^(1-SIGMA)*ni^SIGMA ;


ytot = itot+ctot;


acdb=0;
acdh=0;
ackh=0;
acks=0;
ackn=0;
acloi=0;
aclos=0;
aclon=0;
acloib=0;
aclosb=0;
aclonb=0;

Y_SS = (ys+yn);
ITOT_SS = (itot);
CTOT_SS = (ctot);
Q_SS = (q);
LTOT_SS = loi+lon+los;
NTOT_SS = ntot;
KS_SS = ks;
KN_SS = kn;
KH_SS = kh;
RKH_SS = rkh;
RKS_SS = rks;
RKN_SS = rkn;
D_SS = d;
LOI_SS = loi;
LON_SS = lon;
LOS_SS = los;

intshare = 1;

loitoy=loi/4/(ys+yn);
lostoy=los/4/(ys+yn);
lontoy=lon/4/(ys+yn);
lotoy=(lon+loi+los)/4/(ys+yn);
cbtoy = cb/(ys+yn) ;
cstoy = cs/(ys+yn) ;
cntoy = cn/(ys+yn) ;
chtoy = ch/(ys+yn) ;
citoy = ci/(ys+yn) ;


kab=1-(d/(loi+lon+los));

data_CC=0;
data_KAB=kab;
data_II=0;
data_LLI=0;
data_LLS=0;
data_LLN=0;
data_LLH=0;
data_LO=0;
data_LOI=0;
data_LOS=0;
data_LON=0;
data_QQ=0;
data_YY=0;
data_TFP=0;
data_NN=0;

zks=1;
zkn=1;
zkh=1;

dks=DELTA;
dkn=DELTA;
dkh=DELTA;

CS_SS = cs;
CN_SS = cn;
CE_SS = cs+cn;
RS_SS = rs;
RN_SS = rn;
LS_SS = lambdas;
LN_SS = lambdan;
LB_SS = lambdab;
CB_SS = cb;
CH_SS = ch;

ucb = cb^-SIGMAB;
uci = ci^-SIGMAI;
uch = ch^-SIGMAH;
ucs = cs^-SIGMAS;
ucn = cs^-SIGMAN;

a_bi=0;
a_bs=0;
a_bn=0;
a_bh=0;
a_j=0;
a_k=0;
a_ms=0;
a_mn=0;
a_mi=0;
a_mh=0;
a_p=0;
a_z=0;

e_a_bi=0;
e_a_bs=0;
e_a_bn=0;


xxx = [ ...
a_bi
a_bs
a_bn
a_j
a_k
a_mn
a_ms
a_mi
a_p
a_z
acdb
acdh
acks
ackh
ackn
acloi
aclos
aclon
acloib
aclosb
aclonb
cb
ci
ch
cs
cn
ctot
d
data_TFP
data_YY
data_CC
data_II
data_LO
data_LOI
data_LOS
data_LON
data_QQ
data_LLS
data_LLN
data_NN
dks
dkn
dkh
e_a_bi
e_a_bs
e_a_bn
hi
hh
hs
hn
itot
ks
kn
kh
ktot
lambdab
lambdai
lambdas
lambdan
loi
los
lon
nh
ni
ntot
q
ri
rn
rh
rkh
rkn
rks
rs
rvn
rvs
spread_rn_rh
spread_rs_rh
spread_ri_rh
ucb
uch
uci
ucs
ucn
wh
wi
ys
yn 
zks
zkn
zkh] ;



ys = (xxx) ;

