function F = 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)
% x: VECTOR OF 13 UNKNOWN VARIABLES: c, hn, a, atau, an, pn, ax, am, ab, al, wn, kn, in
% GIVEN PARAMETERS AND KNOWN VARIABLES
% UNKNOWNS: 13 variables
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);
% SYSTEM OF EQUATIONS
F = zeros(13,1);
F(1) = (c - (hm^(OMEGAM)/OMEGAM + hb^(OMEGAB)/OMEGAB + hl^(OMEGAL)/OMEGAL +hn^(OMEGAN)/OMEGAN))^(-SIGMA)*(1-BETA*(1+r)*ptau);
F(2) = (hn)^(OMEGAN-1) - wn;
F(3) = a - (CHI*atau^(1-1/MUNTAU) + (1-CHI)*an^(1-1/MUNTAU))^(1/(1-1/MUNTAU));
F(4) = pn - (1-CHI)*(CHI*atau^(1-1/MUNTAU) + (1-CHI)*an^(1-1/MUNTAU))^(1/(MUNTAU-1))*(an)^(-1/MUNTAU);
F(5) = pm - ptau*(1-CHIX)*(CHIX*ax^(1-1/MUX) + (1-CHIX)*am^(1-1/MUX))^(1/(MUX-1))*(am)^(-1/MUX);
F(6) = ax - (ab)^(CHIB)*(al)^(1-CHIB);
F(7) = pb - px*CHIB*(ab/al)^(CHIB-1);
F(8) = pl - px*(1-CHIB)*(ab/al)^(CHIB); 
F(9) = wn - pn*AN*(1-ALPHAN)*(kn/hn)^(ALPHAN);
F(10) = c + im + ib + il + in - a;
F(11) = an - AN*(kn)^(ALPHAN)*(hn)^(1-ALPHAN);
F(12) = atau - (CHIX*ax^(1-1/MUX) + (1-CHIX)*am^(1-1/MUX))^(1/(1-1/MUX));
F(13) = in - DELTA*kn;
end