function [ns1, nf1, Ns1] = fveq(ns0, nf0, Ns0, rk, ALPHA, BETA, DELTA, EPS, PHI, SIGMA, VARPHI, GAMMA, Zg, Zm, TM)
options = optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

function f = lhs_sol(x)
    ns = x(1); 
    nf = x(2);
    Ns = x(3);

    kpns = (ns*(1 + ns)^GAMMA*TM)/((1 - ns*TM)^2*Zm) + (GAMMA*ns*(1 + ns)^(-1 + GAMMA))/((1 - ns*TM)*Zm) + ...
  (1 + ns)^GAMMA/((1 - ns*TM)*Zm);

    kpnf = (nf*(1 + nf)^GAMMA*TM)/((1 - nf*TM)^2*Zm) + (GAMMA*nf*(1 + nf)^(-1 + GAMMA))/((1 - nf*TM)*Zm) + ...
  (1 + nf)^GAMMA/((1 - nf*TM)*Zm);

    implied_yd_ns = rk * kpns / (SIGMA/ns * PHI/(1+ns) * (EPS/(EPS-1))^-2);

    implied_yd_nf = rk * kpnf / (SIGMA/nf * PHI/(1+nf) * (EPS/(EPS-1))^-2);

% ps is the expoenential ratio between ns nf
    ps = (nf^(1/(-1 + EPS))*(1 + nf)^PHI)/(ns^(1/(-1 + EPS))*(1 + ns)^PHI*...
      ((nf^(1/(-1 + EPS))*(1 + nf)^PHI)/(ns^(1/(-1 + EPS))*(1 + ns)^PHI))^SIGMA);

% ps is also the marked up marginal cost
    w = (((1 - EPS)^2*ns^(1/(-1 + EPS))*(1 + ns)^PHI*ps*Zg^(1 - ALPHA))/...
       (((-1 + 1/BETA + DELTA)/ALPHA)^ALPHA*EPS^2))^(1/(1 - ALPHA)) - ...
     ALPHA*(((1 - EPS)^2*ns^(1/(-1 + EPS))*(1 + ns)^PHI*ps*Zg^(1 - ALPHA))/...
        (((-1 + 1/BETA + DELTA)/ALPHA)^ALPHA*EPS^2))^(1/(1 - ALPHA));

% now we have marginal cost we calculate aggregate production     
    ydprod = ((1 + ns)^PHI*Ns*ps*((ALPHA*w)/((1 - ALPHA)*(-1 + 1/BETA + DELTA)))^ALPHA*Zg^(1 - ALPHA))/SIGMA;

% finally aggregate demand is
    cd = w/((Ns - ((1 + ns)^PHI*Ns*ps*(-1 + SIGMA))/((1 + nf)^PHI*(ps^(-SIGMA))^(1/(1 - SIGMA))*SIGMA))^VARPHI*(2 - BETA));

    Nd = Ns - ((1 + ns)^PHI*Ns*ps*(-1 + SIGMA))/((1 + nf)^PHI*SIGMA);
    kOVn = (w/rk * ALPHA/(1-ALPHA));
    kg = kOVn * Nd;
    km = Zm^-1 * (((1+ns)^GAMMA/1-ns * TM) + ((1+nf)^GAMMA/1-nf * TM));

    k = kg + km;

    yddmnd = cd + DELTA * k;


%implied yd is equal for ns nf
    f(1) = implied_yd_nf - implied_yd_ns;

% implied yd is equal to production yd
    f(2) = implied_yd_ns - ydprod;

% implied yd is equal to aggregate demand
    f(3) = implied_yd_nf - yddmnd;

end

cat = fsolve(@lhs_sol, [ns0, nf0, Ns0], options);

ns1 = cat(1);
nf1 = cat(2);
Ns1 = cat(3);

end