This is an exercise from the book "A First Course in Random Matrix Theory" by Potters and Bouchaud.
We're interested in finding the eigenvalue density $\rho(\lambda)$ of a large random matrix $\mathbf{A}$ whose $k$th normalized moment is given by $\tau(\mathbf{A}^k)=\lim_{N\to \infty}\frac{1}{N}{\rm{tr}}(\mathbf{A}^k)=\frac{1}{k}$ (assuming that in the large $N$ limit all these values are deterministic).
I start by finding the series expansion for Stieltjes transform $g(z) = \int \rho(x)/(z-x) dx$. We know that for large $z$ away from the real axis, $g(z)$ has a series expansion $g(z)=\sum_{k=0}^\infty \frac{1}{z^{k+1}}\tau(\mathbf{A}^k)$. Therefore, for our case $$g(z) = \frac{1}{z} + \sum_{k=1}^\infty \frac{1}{z^{k+1}}\frac{1}{k}.$$ Summing up the series, I find that $g(z)=\frac{1}{z}-\frac{\log \left(1-\frac{1}{z}\right)}{z}.$ Now, to get $\rho(x)$, we use the identity $\rho(x) = \frac{1}{\pi}\lim_{\epsilon\to 0^+} \Im g(x-i\epsilon)$. Doing so, I find that $$\Im g(x-i\epsilon)=-\frac{x \arg \left(1-\frac{1}{x-i \epsilon }\right)}{x^2+\epsilon ^2},$$ which implies that $\rho(x)= 1/x$ for for $0<x<1$ and $0$ otherwise.
But, of course, the problem is that this is not a valid density as it diverges at $0$. If we pretend that it's valid, it reproduces all the moments except the $0$th one correctly.
So, my question is what am I doing wrong here? And what is the correct density $\rho(x)$ for these matrices?
Your computations are correct, and I believe that the problem is incorrect.
Indeed, we can show that there is no finite measure supported on $[0, 1]$ whose $k$th moments are $1/k$. This is related to the Hausdorff moment problem, the resolution of which states that a sequence $m$ is a moment sequence for a Radon measure if it is completely monotonic, ie. that the $k$th order difference sequence satisfies $$ (-1)^k(\Delta^k m)_n \geq 0 $$ for all $k, n \geq 0$.
Now consider a sequence $m = (m_0, 1, 1/2, 1/3, \dotsc, )$ for some $m_0 \in (0, \infty)$. The zeroth term of the corresponding difference sequence is \begin{align*} (\Delta^k m)_0 = \sum_{j=0}^k (-1)^j \binom{k}{j} m_{k-j} = (-1)^k \Bigl[\sum_{j = 1}^k (-1)^j \binom{k}{j} \frac{1}{j} + m_0\Bigr]. \end{align*}
But there is a very nice identity that states that $$ \sum_{j = 1}^k (-1)^j \binom{k}{j} \frac{1}{j} = -H_k, $$ where $H_k = \sum_{n=1}^k 1/n$ is the harmonic series. Thus, for $m$ to be completely monotone, we require that, for all $k$, $$ (-1)^k(\Delta^k m)_0 = -H_k + m_0 \geq 0. $$
But since $H_k$ diverges, this is impossible. That is, there is no finite $m_0$ that makes $m$ a difference sequence, and so no finite measure $\mu$ on $[0, 1]$ can have moments $\mu(x^k) = 1/k$.