ssigma=2;
ppsi=0.5;
oomegac=0.5;
tthetac=0.5;
ggamma=0.5;
pphi=3;
oomegal=0.5;
tthetal=1.1;
xxil=1;
Ah=1;
aalphah=0.7;
An=1;
aalphan=0.7;
xxi=0.1;
bbeta=0.98;
xxib=0.01;
Bs=0.5;
nnu=-log(bbeta);
zzetais = 0.7;
zzetaphs = 0.7;
zzetapfs = 0.7;
zzetarho = 0.7;
xxil = 2;
pphiw = 1;
pphip = 1;
eepsilonw = 2;
eepsilonn = 2;

x0 = zeros(13,1);
x = fsolve(@(x)steadystate(x,ssigma,ppsi,oomegac,tthetac,ggamma,pphi,oomegal,tthetal,xxil,Ah,aalphah,An,aalphan),x0);

Css  = exp(x(1));
CTss = exp(x(2));
CNss = exp(x(3));
CHss = exp(x(4));
CFss = exp(x(5));
Lss  = exp(x(6));
LHss = exp(x(7));
LNss = exp(x(8));
PTss = exp(x(9));
PNss = exp(x(10));
Ess  = exp(x(11));
WHss = exp(x(12));
WNss = exp(x(13));

YHss = Ah*LHss^aalphah;

CT = oomegac^(1/tthetac)*CTss^((tthetac-1)/tthetac)/Css^((tthetac-1)/tthetac);
CN = (1-oomegac)^(1/tthetac)*CNss^((tthetac-1)/tthetac)/Css^((tthetac-1)/tthetac);
CH = ggamma^(1/pphi)*CHss^((pphi-1)/pphi)/CTss^((pphi-1)/pphi);
CF = (1-ggamma)^(1/pphi)*CFss^((pphi-1)/pphi)/CTss^((pphi-1)/pphi);
ELLH = oomegal^(1/tthetal)*LHss^((tthetal-1)/tthetal)/Lss^((tthetal-1)/tthetal);
ELLN = (1-oomegal)^(1/tthetal)*LNss^((tthetal-1)/tthetal)/Lss^((tthetal-1)/tthetal);
uf = Css^(-ssigma+1/tthetac)*oomegac^(1/tthetac)*CTss^(1/pphi-1/tthetac)*(1-ggamma)^(1/pphi)*CFss^(-1/pphi);

lmat = [1 0 0 0 eepsilonw*aalphah/pphiw*(1/tthetac-ssigma) eepsilonw*aalphan/pphiw*(1/tthetac-ssigma) 0 0 0 0;
    -CT 1 0 0 eepsilonw*aalphah/pphiw*(1/pphi-1/tthetac) 0 0 0 0 0;
    -CN 0 0 0 0 -eepsilonw*aalphan/pphiw*1/tthetac 0 0 0 1;
    0 -CH 0 1/pphi -eepsilonw*aalphah/pphiw*1/pphi 0 0 0 CHss/CFss 0;
    0 -CF 0 -1/pphi 0 0 0 0 1 0;
    0 0 1 0 -eepsilonw*aalphah/pphiw*(ppsi+1/tthetal) -eepsilonw*aalphan/pphiw*(ppsi+1/tthetal) 0 0 0 0;
    0 0 -ELLH 0 eepsilonw*aalphah/pphiw*1/tthetal 0 (aalphah-1) 0 -YHss/CFss*aalphah 0;
    0 0 -ELLN 0 0 eepsilonw*aalphan/pphiw*1/tthetal 0 (eepsilonn-1)/pphip*(aalphan-1) 0 -aalphan;
    0 0 0 0 0 eepsilonw*aalphan/pphiw 0 -(eepsilonn-1)/pphip 0 0;
    0 0 0 0 eepsilonw*aalphah/pphiw 0 -1 0 0 0;];

consmat = [Css^(1-ssigma);0;0;0;0;-Lss^(1+ppsi);0;0;0;0;];

lambdas = lmat^(-1)*consmat;

l1ss = lambdas(1);
l2ss = lambdas(2);
l3ss = lambdas(3);
l7ss = lambdas(4);
l10ss = lambdas(5);
l11ss = lambdas(6);
l12ss = lambdas(7);
l13ss = lambdas(8);
l14ss = lambdas(9);
l15ss = lambdas(10);

l4ss = 0;
l5ss = 0;
l6ss = 0;
l8ss = 0;
l9ss = 0;
l16ss = 0;
l17ss = 0;
l18ss = 0;