I can get,
$$h_n \approx \ln n+\gamma$$
By, saying
$$\frac{h_{n}-h_{n-1}}{n-(n-1)}=\frac{1}{n}$$
Then letting $h_n=h(n)$ be continuous and differentiable on $(0,\infty)$. So that,
$$h'(n) \approx \frac{1}{n}$$
And thus,
$$h_n \approx \ln n+C$$
Where we will let $C$ give us the best approximation as $n \to \infty$ that is,
$$C=\lim_{n \to \infty} (h_n-\ln n)$$
What about the other terms, I'm curious to see how they are derived.
Let $x_n=H_n-(\ln n+\gamma)$ with $H_n=\sum_{k=1}^n\frac{1}{k}$ and evaluate the following limit by using the Stolz-Cesaro Theorem ($\ast/0$ case), $$\lim_{n\to \infty}\frac{x_n}{1/n}=\lim_{n\to \infty}\frac{x_{n+1}-x_n}{1/(n+1)-1/n}=\lim_{n\to \infty}\frac{1/(n+1)-\ln(1+1/n)}{1/(n+1)-1/n}=\frac{1}{2}.$$ Hence $H_n=\ln(n)+\gamma+\frac{1}{2n}+o(1/n)$.
More terms can be found in a similar way (a general formula involves the Bernoulli numbers). Another method to find the asymptotic expansion of the harmonic numbers employs Euler–Maclaurin formula.
$\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} \lim_{n \to \infty}\bracks{H_{n} - \ln\pars{n}} & = \lim_{n \to \infty}\braces{\sum_{k = 1}^{n}{1 \over k} + \bracks{\ln\pars{1 \over 2} + \ln\pars{2 \over 3} + \cdots + \ln\pars{n - 1 \over n}}} \\[5mm] & = \lim_{n \to \infty}\braces{1 + \sum_{k = 2}^{n}\bracks{{1 \over k} + \ln\pars{k - 1 \over k}}} \\[5mm] & = \lim_{n \to \infty}\bracks{1 + \sum_{k = 2}^{n}\int_{0}^{1}\pars{{1 \over k} + {1 \over t - k}}\,\dd t} \\[5mm] & = 1 + \int_{0}^{1}\sum_{k = 0}^{\infty} \pars{{1 \over k + 2} - {1 \over k + 2 - t}}\,\dd t = 1 + \int_{0}^{1}\pars{H_{1 - t} - H_{1}}\,\dd t \\[5mm] & = \int_{0}^{1}H_{1 - t}\,\dd t = \int_{0}^{1}\bracks{-\,\totald{\ln\pars{\Gamma\pars{2 - t}}}{t} + \gamma}\,\dd t = \bbx{\gamma} \end{align}
