function [c,hn,kn,a,atau,an,pn,ax,am,ab,al,wn,in] = mxn_steadystate_helper(BETA,DELTA,SIGMA,OMEGAB,OMEGAL,OMEGAM,OMEGAN,CHIB,ALPHAN,AN,CHI,CHIX,MUX,MUNTAU,pb,pl,px,pm,r,hm,hb,hl,im,ib,il,ptau)
KM_O_HM = 11.5315;
c0      = 150;
hn0     = 11.1697;
kn0     = hn0 * KM_O_HM;
atau0   = 30;
an0     = AN * (hn0)^(ALPHAN) * (kn0)^(1-ALPHAN);
pn0     = (1-CHI)/CHI * (an0/atau0)^(-1/MUNTAU);
ax0     = (px/(ptau*CHIX*atau0^(-1/MUX)))^(-MUX);
am0     = (pm/(ptau*(1-CHIX)*atau0^(-1/MUX)))^(-MUX);
ab0     = 30;
al0     = ab0 * (pb/(px*CHIB))^(-1/(CHIB-1));
wn0     = (hn0)^(OMEGAN-1);
in0     = DELTA * kn0;
%a0      = c0 + in0 + im + il + ib;
a0      = (CHI*atau0^(1-1/MUNTAU) + (1-CHI)*an0^(1-1/MUNTAU))^(1/(1-1/MUNTAU));

x0 = [c0,hn0,kn0,a0,atau0,an0,pn0,ax0,am0,ab0,al0,wn0,in0];
F = @(x) myfun(x,BETA,DELTA,SIGMA,OMEGAB,OMEGAL,OMEGAM,OMEGAN,CHIB,ALPHAN,AN,CHI,CHIX,MUX,MUNTAU,pb,pl,px,pm,r,hm,hb,hl,im,ib,il,ptau);
max = 1000000;
options = optimoptions(@fsolve, 'MaxFunctionEvaluations', max, 'MaxIterations', max);
x = fsolve(F,x0,options);
%OUTPUTS
c    = x(1);
hn   = x(2);
kn   = x(3);
a    = x(4);
atau = x(5);
an   = x(6);
pn   = x(7);
ax   = x(8);
am   = x(9);
ab   = x(10);
al   = x(11);
wn   = x(12);
in   = x(13);

end