  It is known that the generating function $f$ of a sequence of Toeplitz
matrices $\{T_n(f)\}_n$ may not describe the asymptotic distribution of the
eigenvalues of $T_n(f)$ if $f$ is not real. In this paper, we assume as a
working hypothesis that, if the eigenvalues of $T_n(f)$ are real for all $n$,
then they admit an asymptotic expansion of the same type as considered in
previous works [1,10,12,13], where the first function $g$ appearing in this
expansion is real and describes the asymptotic distribution of the eigenvalues
of $T_n(f)$. After validating this working hypothesis through a number of
numerical experiments, drawing inspiration from [12], we propose a matrix-less
algorithm in order to approximate the eigenvalue distribution function $g$. The
proposed algorithm is tested on a wide range of numerical examples; in some
cases, we are even able to find the analytical expression of $g$. Future
research directions are outlined at the end of the paper.
