function h = h_helper(h0,alph,beta,mu)
h = fsolve(@(h) h-alph/(alph + 1 + 1/h^alph*beta - 1/mu), h0);
end
 