%% Solving for the deterministic steady state

function F = steadystate(x,ssigma,ppsi,oomegac,tthetac,ggamma,pphi,oomegal,tthetal,xxil,AH,aalphah,AN,aalphan)

%      1 2  3  4  5  6 7  8  9  10 11 12 13
% x = [C,CT,CN,CH,CF,L,LH,LN,PT,PN,E,WH,WN]

F(1) = exp(x(1)) - (oomegac^(1/tthetac)*exp(x(2))^((tthetac-1)/tthetac) + (1-oomegac)^(1/tthetac)*exp(x(3))^((tthetac-1)/tthetac))^(tthetac/(tthetac-1));
F(2) = exp(x(6)) - (oomegal^(1/tthetal)*exp(x(7))^((tthetal-1)/tthetal) + (1-oomegal)^(1/tthetal)*exp(x(8))^((tthetal-1)/tthetal))^(tthetal/(tthetal-1));
F(3) = exp(x(2)) - (ggamma^(1/pphi)*exp(x(4))^((pphi-1)/pphi) + (1-ggamma)^(1/pphi)*exp(x(5))^((pphi-1)/pphi))^(pphi/(pphi-1));
F(4) = (1-oomegac)^(1/tthetac)*exp(x(3))^(-1/tthetac)/exp(x(10)) - oomegac^(1/tthetac)*exp(x(2))^(-1/tthetac)/exp(x(9));
F(5) = ggamma^(1/pphi)*exp(x(4))^(-1/tthetac) - (1-ggamma)^(1/pphi)*exp(x(5))^(-1/pphi);
F(6) = exp(x(1))^(1/tthetac)*(1-oomegac)^(1/tthetac)*exp(x(3))^(-1/tthetac)/exp(x(10)) -1;
F(7) = (1-oomegac)^(1/tthetac)*exp(x(3))^(-1/tthetac)/exp(x(10)) - oomegac^(1/tthetac)*exp(x(2))^(1/pphi-1/tthetac)*(1-ggamma)^(1/pphi)*exp(x(5))^(-1/pphi)/exp(x(11));
F(8) = xxil*exp(x(6))^(ppsi+1/tthetal)*oomegal^(1/tthetal)*exp(x(7))^(-1/tthetal)/(exp(x(1))^(1/tthetac-ssigma)*oomegac^(1/tthetac)*exp(x(2))^(1/pphi-1/tthetac)*ggamma^(1/pphi)*exp(x(4))^(-1/pphi)) - exp(x(12))/exp(x(11));
F(9) = xxil*exp(x(6))^(ppsi+1/tthetal)*(1-oomegal)^(1/tthetal)*exp(x(8))^(-1/tthetal)/(exp(x(1))^(1/tthetac-ssigma)*(1-oomegac)^(1/tthetac)*exp(x(3))^(-1/tthetac)) - exp(x(13))/exp(x(11));
F(10) = aalphah*AH*exp(x(7))^(aalphah-1) - exp(x(12))/exp(x(11));
F(11) = aalphan*AN*exp(x(8))^(aalphan-1) - exp(x(13))/exp(x(10));
F(12) = exp(x(5)) + exp(x(4)) - AH*exp(x(7))^aalphah;
F(13) = exp(x(3)) - AN*exp(x(8))^aalphan;