function [err, q, mu, Omega, Z, W, Lp, K, B, U, Uq, C, Mu, P]  = findequilibrium(X, w, z, p, type, steady, Nfirms)

D       = X(1); 

if strcmp(steady, 'old')

Y       = 1; 
N       = 1; 

else
  
Y       = X(2); 
N       = X(3); 
    
end

if nargin == 7
   
    N = Nfirms;     % this is for computing aggregate productivity and markups as a function of number of firms
    
end

R       = p.r + p.delta;

if strcmp(type, 'market')
 
    mu           = ones(size(z)); 
    
    q            = zeros(size(z));

    good         = z >= p.sigma/(p.sigma - 1)*exp(-1/p.epsi)*D; 

    mu(good)     = bisect('firmfoc', 1, 250000, z(good), D, p); 

    q(good)      = p.sigma^(p.sigma/p.epsi)*(1 - 1./mu(good)).^(p.sigma/p.epsi); 
    
    mu           = mu./(1 + p.xis);      % this is markup over true marginal cost (excl. subsidy)
    
else
    
    q            = (1 - min(p.epsi*log(p.sigma/(p.sigma - 1)*D*(1 + p.taus)./z), 1)).^(p.sigma/p.epsi); 

    mu           = ones(size(z))*(1 + p.taus); 
    
end
    
U       = upsilon(q, p); 
Uq      = (p.sigma - 1)/p.sigma*exp((1 - q.^(p.epsi/p.sigma))/p.epsi);                 % U'(q)

Z       = (N*w'*(q./z))^(-1);

Mu      = (w'*(1./mu.*Uq.*q)/(w'*(Uq.*q))).^(-1);                                 % aggregate markup


P       = (N*w'*(Uq.*q))^(-1);       % Demand index so we can compute prices


Omega   = Z/Mu; 

W       = (1 - p.alpha)*(((Omega^(1 - p.theta) - (1 - p.phi))/p.phi)^(1/(1 - p.theta))*(R/p.alpha)^(-p.alpha))^(1/(1 - p.alpha));        % Invert the CES Cost aggregator                                              % Labor only factor of production                

Lp      = (1 - p.alpha)*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/W; 
K       =       p.alpha*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/R;
B       = (1 - p.phi)/Mu*(1/Omega)^(1 - p.theta)*Y;

C       = Y - p.delta*K - B;

err     = zeros(numel(X), 1); 

err(1)  = N*w'*U - 1;

if strcmp(steady, 'new')

    err(2)    = W - p.psi*C*(Lp + p.F*p.varphi*N)^p.nu; 

    
    if strcmp(type, 'market')
   
    err(3)    = p.F*W - (1 + p.xie)/(p.beta^(-1) - 1 + p.varphi)*Y/N*(1 + p.xis - 1/Mu);  

    else
        
    err(3)    = p.F*W - 1/(p.beta^(-1) - 1 + p.varphi)*Y/N*(Z/D/(1 + p.taus) - 1)/Mu;      
    
    end

end

