function [ ch,R,ks,ds,div_s,ns,phis,spreads,kc,dc,div_c,nc,phic,spreadc,q,ws,wc,lamda]=solve_fore(BETAH,Z_ss,CHI,ALPHAC,THETA,SIGMAS,SIGMAC,ETA_S,EPSILON,Y_ss,ETA_C,OMEGAK)

R=1/BETAH;

lnr = fsolve(@(x) solvessd(exp(x(1)),x(2),BETAH,Z_ss,CHI,ALPHAC,THETA,SIGMAS,SIGMAC,ETA_S,Y_ss,ETA_C,OMEGAK),[0 0.01], optimset('TolFun',1e-20));
q=exp(lnr(1));
ks=lnr(2);

if ks < 0
lnr = fsolve(@(x) solvessd(exp(x(1)),exp(x(2)),BETAH,Z_ss,CHI,ALPHAC,THETA,SIGMAS,SIGMAC,ETA_S,Y_ss,ETA_C,OMEGAK),[0 log(0.9)], optimset('TolFun',1e-20));
q=exp(lnr(1));
ks=exp(lnr(2));
end

if ks > 1
lnr = fsolve(@(x) solvessd(exp(x(1)),1-exp(x(2)),BETAH,Z_ss,CHI,ALPHAC,THETA,SIGMAS,SIGMAC,ETA_S,Y_ss,ETA_C,OMEGAK),[0 log(1-0.5)], optimset('TolFun',1e-20));
q=exp(lnr(1));
ks=1-exp(lnr(2));
end


 kc=1-ks;
 Rs=(Z_ss+q)/q;
spreads=Rs-R;
phis=R/(THETA-spreads);
ns=ks*q/phis;
ds=q*ks-ns;
div_s=(1-SIGMAS)*(ks*(Z_ss+q)-R*ds);

phic=((Z_ss+q)/R-q-ALPHAC*kc)/(OMEGAK*q)+1/CHI;
nc=kc*q/phic;
wc=Y_ss*ETA_C;
dc=(SIGMAC*(Z_ss+q)*kc - nc +(1-SIGMAC)* wc)/SIGMAC/R;
Rc=(Z_ss+q)/(q+ALPHAC*kc);
spreadc=Rc-R;
div_c=(1-SIGMAC)*(kc*(Z_ss+q)-R*dc);
g=EPSILON*Y_ss;
ch=Y_ss-g-ALPHAC/2*kc^2;
ws=ETA_S*Y_ss;
lamda=R;

 