function [Mu, Z, Eprof] = solution(z, p) 


[mu, ~, ~, ~, ~, ~, Mu, Z]  = industryequilibrium([], z, p);

nz          = size(p.z, 1); 
S           = size(z,   2); 

Profz       = zeros(nz, 1); 

muold       = [repmat(p.gamma/(p.gamma - 1), 1, S); mu]; 


parfor i = 1  : nz
    
   znew          = [repmat(p.z(i), 1, S); z]; 
    
   [~, profnew]  = industryequilibrium(muold, znew, p);
    
   Profz(i)      = mean(profnew(1, :)); 
   
end

Eprof            = (p.w'*Profz);

