  Let $ S $ be a positivity preserving symmetric linear operator acting on
bounded functions. The nonlinear equation $ -\frac{1}{m}=z+Sm $ with a
parameter $ z $ in the complex upper half-plane $ \mathbb{H} $ has a unique
solution $ m $ with values in $ \mathbb{H} $. We show that the $ z $-dependence
of this solution can be represented as the Stieltjes transforms of a family of
probability measures $ v $ on $ \mathbb{R} $. Under suitable conditions on $ S
$, we show that $ v $ has a real analytic density apart from finitely many
algebraic singularities of degree at most three.
  Our motivation comes from large random matrices. The solution $ m $
determines the density of eigenvalues of two prominent matrix ensembles; (i)
matrices with centered independent entries whose variances are given by $ S $
and (ii) matrices with correlated entries with a translation invariant
correlation structure. Our analysis shows that the limiting eigenvalue density
has only square root singularities or a cubic root cusps; no other
singularities occur.
