function [mu, prof, omega, Mus, Zs, Omegas, Mu, Z, Omega]  = industryequilibrium(muold, z, p)
  
if isempty(muold)

    mu      = vec(p.gamma/(p.gamma - 1)*ones(size(z)));      % guess all charge lowest markup
    
else
    
    mu      = muold(:); 
    
end

clear equilibrium

for it = 1 : 5
   
    mu = equilibrium(mu, z, p); 
   
end

% Anderson

ftarget          = @(mu) equilibrium(mu, z, p); 

parand.theta     = 0.01;
parand.tau       = 0.001;
parand.D         = 1e6;
parand.epsilon   = 1e-6;
parand.mem_size  = 5;
parand.itermax   = 150; 
parand.tol       = 1e-4;

mu               = andersen(mu(:), ftarget, parand); 

mu               = reshape(mu, size(z, 1), size(z, 2)); 

omega            = (mu./z).^(1 - p.gamma)./sum((mu./z).^(1 - p.gamma));  
    
% Sectoral variables

Mus              = sum(omega./mu).^(-1);                                                    % sectoral markup

Zs               = sum((mu./Mus).^(-p.gamma).*z.^(p.gamma - 1)).^(1/(p.gamma - 1)); 

Omegas           = (Mus./Zs).^(1 - p.sigma)./mean((Mus./Zs).^(1 - p.sigma)); 

Mu               = mean(Omegas./Mus).^(-1);                                                 % aggregate markup

Z                = mean((Mus./Mu).^(-p.sigma).*Zs.^(p.sigma - 1)).^(1/(p.sigma - 1));

Omega            = Z/Mu; 


prof             = (1 - 1./mu).*omega.*Omegas;  % careful: this is not scaled by Y or the subsidy
