function [mu, err] = equilibrium(muold, z, p)

persistent it


if isempty(it)
      
     it = 1;
        
end


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

omega   = (mu./z).^(1 - p.gamma);           % Faster than 

omega   = omega./sum(omega);                % omega = (mu./z).^(1 - p.gamma)./sum((mu./z).^(1 - p.gamma));  
 
    
epsi    = max((omega./p.sigma + (1 - omega)./p.gamma).^(-1), 1.1); 
    
mu      = epsi./(epsi - 1);


mu      = mu(:); 
    
err     = max(abs((mu - muold)./muold));

%fprintf('%4i %6.2e \n', [it, err]);       % report difference

mu      = muold + 1/5*(mu - muold);        % update guess

it      = it + 1;
