function [ns1] = nssol(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(ns)

% Caluculate implied Ns form ns
    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);

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

% Calculate wages that set ps = 1
    mc = (EPS/(EPS-1))^-2 * ns^(1/(EPS-1)) * (1+ns)^PHI;
    w =-((-1 + ALPHA)*((mc*Zg^(1 - ALPHA))/(rk/ALPHA)^ALPHA)^(1/(1 - ALPHA)));
    Ns = implied_yd_ns/(2*ns^(EPS/(EPS - 1))*(1 + ns)^PHI*Zg^(1 - ALPHA)*((w/rk)*(ALPHA/(1 - ALPHA)))^ALPHA);
  

% Calculate output production
    ydprod = 2*Ns*ns^(EPS/(EPS - 1))*(1 + ns)^PHI*Zg^(1 - ALPHA)*((w/rk)*(ALPHA/(1 - ALPHA)))^ALPHA; 
    
% Calculate output demand
    Nd = 2 * Ns;
    cd = w/((2 - BETA)*Nd^VARPHI);
    cd = w/((2*Ns)^VARPHI*(2 - BETA)); 
    kg = ((w/rk)*(ALPHA/(1 - ALPHA)))*2*Ns; 
    km = (1 + ns)^GAMMA/((1 - ns*TM)*Zm) + (1 + ns)^GAMMA/((1 - ns*TM)*Zm); 
    k = kg+km;
    inv = DELTA * k; 
    yddmnd = cd + inv; 

f = yddmnd - ydprod;
end

cat = fsolve(@lhs_sol, ns0, options);

ns1 = cat(1);

end