A better approximation of $H_n $ I'm convinced that
$$H_n \approx\log(n+\gamma) +\gamma$$ is a better approximation of the $n$-th harmonic number than the classical $$H_n \approx \log(n) +\gamma$$ 
Specially for small values of $n$. I leave some values and the error:


Just to make things clearer, I calculate the value between two numbers as follows.
Say $n$ is the optimal and $a$ is the apporximation, then $E = \frac{n-a}{n}$. $L_1$ stands for my approximation and $L_2$ for the classical one, and the errors $E_2$ and $E_1$ correspond to each of those (I mixed up the numbers).
It is clear that this gives an over estimate but tends to the real value for larger $n$.
So, is there a way to prove that the approximation is better?

NOTE: I tried using the \begin{tabular} environment but nothing seemed to work. Any links on table making in this site?

The asymptotic expansion of the Harmonic numbers $H_n$ is given by
$$\log n+\gamma+\frac{1}{2n}+\mathcal{O}\left(\frac{1}{n^2}\right).$$
The Maclaurin series expansion of the natural logarithm tells us $\log(1+x)=x+\mathcal{O}(x^2)$, and we can use this in your formula by writing $\log(n+\epsilon)=\log n+\log(1+\epsilon/n)$ and expanding:
$$\log(n+\gamma)+\gamma=\log n+\gamma\;\;\;+\frac{\gamma}{n}+\mathcal{O}\left(\frac{1}{n^2}\right).$$
Your approximation is asymptotically better than the generic one because the $\gamma=0.577\dots$ in its expansion is closer to the true coefficient $\frac{1}{2}$ than the illicit coefficient $0$ in the generic formula given by the usual $H_n\sim \log n +\gamma+0/n$. This also explains why it is asymptotically an over estimation.

As marty said in his answer, the expansion comes from the Euler-Maclaurin formula:
$$\sum_{n=a}^b f(n)=\int_a^b f(x)dx+\frac{f(a)+f(b)}{2}+\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right).$$
Here we let $a=1,b=n$ (rewrite the index to a different letter) and $f(x)=1/x$.
