function F = NKgertlerfun_replfix_grossinvdo(unknown,param)

gs=param(1);
beta=param(2);
sigma=param(3);
rho=param(4);
xi=param(6);
bs=param(7);
eta=param(8);
delta=param(9);
alpha=param(10);
v1=param(11);
v2=param(12);
v3=param(13);
omega=param(14);
nw =param(15);
x=param(16);
gamma=param(17);
theta=param(18);
psi=param(19);
es=param(22);


bigomega =    unknown(1);
epsilon  =    unknown(2);
cw       =    unknown(3);
cr       =    unknown(4);
pi       =    unknown(5);
lamda    =    unknown(6);
hr       =    unknown(7);
hw       =    unknown(8);
lw       =    unknown(9);
lr       =    unknown(10);
dr       =    unknown(11);
dw       =    unknown(12);
c        =    unknown(13);
m        =    unknown(14);
a        =    unknown(15);
t        =    unknown(16);
l        =    unknown(17);
r        =    unknown(18);
i        =    unknown(19);
inv      =    unknown(20);
k        =    unknown(21);
rk       =    unknown(22);
w        =    unknown(23);
e        =    unknown(24);

%bigomega epsilon cw cr pi lamda hr hw lw lr dr dw c m a t l r i inv k rk   

mc=(theta-1)/theta;
y=(k/((1+nw)*(1+x)))^(1-alpha)*l^alpha;
mu=e/w*psi;

F = [bigomega-omega-(1-omega)*(epsilon^(1/(1-sigma)))*(1/xi)^v3;
    
     epsilon*pi-1+((1/(1+x))^(v3*rho*sigma))*(beta^sigma)*((1+r)^(sigma-1))*gamma;
     pi-1+((1/(1+x))^(v3*rho*sigma))*(beta^sigma)*(((1+r)*bigomega)^(sigma-1));
     
     cr*(1+(v2/v1))-epsilon*pi*((1+r)*lamda*a/((1+nw)*(1+x))+hr);
     cw*(1+v2/v1) - pi* ( ((1+r)*(1-lamda)*a /((1+nw)*(1+x))) +hw );
     c-cw-cr*psi;
     
     lamda*a*(1-omega*(1-epsilon*pi)*(1+r)/((1+nw)*(1+x))) - (1-omega)*a - omega*psi*(dr -epsilon*pi*hr);
     
     lw-1+(v3/v1)*(cw/w);
     lr-1+(v3/v1)*(cr/(xi*w));
     l-lw-lr*psi*xi;
     
     m-((1+i)/i)*(v2/v1)*c;
     a-(m/(1+i))-bs*y-k;
     
     (t-e)*((1+nw)*(1+x))-bs*y*((1+r)-(1+nw)*(1+x))+(1+nw)*(1+x)*(m-gs*y)-m;
     gs*y-y+c+inv;
     
     dr-w*lr*xi-e/psi;
     dw-w*lw-(1/theta)*y+t;
     hr*(1-gamma*((1+x)/(1+r)))-dr;   
     hw*(1-(omega/bigomega)*((1+x)/(1+r)))-dw-(1-omega/bigomega)*((1+x)/(1+r))*hr; %%my version
     %hw*(1-(omega/bigomega)*((1+x)/(1+r)))-dw-(1-omega/bigomega)*((1+x)/(1+r))*hr*psi; %corresponds to old version
     %hw*(1-(omega/bigomega)*((1+x)/(1+r)))-dw-(1/bigomega)*(1-omega)*(epsilon^(1/(1-sigma)))*(1/xi)^v3*((1+x)/(1+r))*hr*psi; %original

     rk-r-delta;
     r-i; 
     inv-(1-(1-delta)/((1+nw)*(1+x)))*k
     
     e/y-es; %???????????????
     mc-(w/alpha)^alpha*(rk/(1-alpha))^(1-alpha);
     k/((1+nw)*(1+x)*l)-w*(1-alpha)/(alpha*rk)];