Summary: We introduce a method to numerically compute equilibrium measures for problems with attractive-repulsive power law kernels of the form \(K(x-y) = \frac{|x-y|^\alpha}{\alpha}-\frac{|x-y|^\beta}{\beta}\) using recursively generated banded and approximately banded operators acting on expansions in ultraspherical polynomial bases. The proposed method reduces what is naïvely a difficult to approach optimization problem over a measure space to a straightforward optimization problem over one or two variables fixing the support of the equilibrium measure. The structure and rapid convergence properties of the obtained operators result in high computational efficiency in the individual optimization steps. We discuss stability and convergence of the method under a Tikhonov regularization and use an implementation to showcase comparisons with analytically known solutions as well as discrete particle simulations. Finally, we numerically explore open questions with respect to existence and uniqueness of equilibrium measures as well as gap forming behaviour in parameter ranges of interest for power law kernels, where the support of the equilibrium measure splits into two intervals.