How to make approximation of rotational partition function of diatomic linear molecules?
Sorry for the late answer. I just discovered this question.

Preludium

Firstly, I think you might have an error or some non-standard notation in your formula for $Q_{\mathrm{rot}}$. The rotational energy levels are given by

\begin{align}
  E_{\mathrm{rot}} = \frac{\hbar^2}{2 I} J ( J + 1 ) \ ,
\end{align}

where $I$ is the moment of inertia.
Then the rotational constant $B_{0} = \frac{\hbar^2}{2 I}$ (in units of energy) is introduced, so

\begin{align}
  E_{\mathrm{rot}} = B_{0} J ( J + 1 ) \ .
\end{align}

Thus

\begin{align}
  Q_{\mathrm{rot}} = \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{B_0 J(J+1)}{kT}\right)} \ ,
\end{align}

whereas your exponent contains another $h$ which should not be there. At least that's the standard notation I am used to. For convenience I'll also introduce the rotational temperature $\Theta_{\mathrm{rot}} = \frac{B_{0}}{k}$, so that

\begin{align}
  Q_{\mathrm{rot}} = \sum_{J=0}^{\infty}{(2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right)} \ .
\end{align}

Main Answer

Ok, now that that's sorted out, let's start the real work. It is actually a pretty standard approximation that is used here: They use the Euler-Maclaurin formula

\begin{align}
  \sum_{n = a}^{b} f(n) &amp;= \int_{a}^{b} f(n) \, \mathrm{d} n + \frac{1}{2} \bigl( f(b) + f(a) \bigr) + \sum_{i=1}^{\infty} (-1)^{i} \frac{b_{2i}}{(2i)!} \Bigl( f^{(2i - 1)}(a) - f^{(2i - 1)}(b) \Bigr)
\end{align}

where $f^{(k)}(a)$ is the $k$th derivative of $f$ evaluated at $a$. The $b_{2i}$s are the absolute values of the even Bernoulli numbers, $b_{2} = \frac{1}{6}$, $b_{4} = \frac{1}{30}$, $b_{6} = \frac{1}{42}$, $...$.
It provides a way to approximate the summation of a function over a discrete variable by an analogous integration over a continuous variable.

In the special case at hand, 

\begin{align}
  f(J) = (2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right) \ ,
\end{align}

where the summation limits are $a=0$ and $b=\infty$ and where $f(\infty) = f^{\prime}(\infty) = \ldots = 0$ this formula boils down to

\begin{align}
  Q_{\mathrm{rot}} = \sum_{J = 0}^{\infty} f(J) &amp;= \int_{0}^{\infty} f(J) \, \mathrm{d} J + \frac{1}{2} f(0) + \sum_{i=1}^{\infty} (-1)^{i} \frac{b_{2i}}{(2i)!} f^{(2i - 1)}(0)
\end{align}

The first term is what is usually called the classical or the high-temperature limit. It can be easily calculated by making the substitution $K = J(J + 1)$ $\Rightarrow$ $\mathrm{d}J = \frac{\mathrm{d}K}{2J + 1}$. So

\begin{align}
  \int_{0}^{\infty} f(J) \, \mathrm{d} J &amp;= \int_{0}^{\infty} (2J+1)\exp\left(-\frac{\Theta_{\mathrm{rot}} J(J+1)}{T}\right) \, \mathrm{d} J \\
  &amp;= \int_{0}^{\infty} \exp\left(-\frac{\Theta_{\mathrm{rot}} K}{T}\right) \, \mathrm{d} K \\
  &amp;= \frac{T}{\Theta_{\mathrm{rot}}}
\end{align}

Now, only the derivatives $f^{(k)}(0)$ need to be sorted out. Those that contain terms below the third order in $\frac{\Theta_{\mathrm{rot}}}{T}$ can be calculated to be:

\begin{align}
  f(0) &amp;= 1 \\
  f^{(1)}(0) &amp;= 2 - \frac{\Theta_{\mathrm{rot}}}{T} \\
  f^{(3)}(0) &amp;= - 12 \frac{\Theta_{\mathrm{rot}}}{T} + 12 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \\
  f^{(5)}(0) &amp;= 120 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - 180 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} + 30 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{4} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{5} \ .
\end{align}

Doing this by hand is very tedious. Best you use something like Mathematica. Plugging all this into the Euler-MacLaurin formula and dropping all terms above second order in $\frac{\Theta_{\mathrm{rot}}}{T}$ you get exactly the result of Gordy &amp; Cook

\begin{align}
  Q_{\mathrm{rot}} &amp;\approx \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{2} \cdot 1 + (-1) \frac{1}{6 \cdot 2!} \left( 2 - \frac{\Theta_{\mathrm{rot}}}{T} \right) \\
 &amp;\quad + (-1)^{2} \frac{1}{30 \cdot 4!} \left( - 12 \frac{\Theta_{\mathrm{rot}}}{T} + 12 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \\
  &amp;\quad + (-1)^{3} \frac{1}{42 \cdot 6!} \left(120 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - 180 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} + 30 \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{4} - \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{5} \right) \\
  &amp;= \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{2} - \frac{1}{6} + \frac{1}{12} \frac{\Theta_{\mathrm{rot}}}{T} - \frac{1}{60} \frac{\Theta_{\mathrm{rot}}}{T} \\
  &amp;\quad + \frac{1}{60} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} - \frac{1}{252} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \\
  &amp;= \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{3} + \frac{1}{15} \frac{\Theta_{\mathrm{rot}}}{T} + \frac{4}{315} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right)
\end{align}

As a side note: If you have a homonuclear diatomic molecule you need to introduce the so called symmetry number $\sigma$ which represents the number of indistinguishable orientations a molecule can have. It originates from symmetry constrains: for homonuclear diatomic molecules there is a two-fold axis of symmetry perpendicular to the internuclear axis which gives rise to two indistinguishable orientations, so the above formula for $Q_{\mathrm{rot}}$ overcounts the number of available quantum states
by a factor of two because of the inherent indistinguishability of the nuclear pair. Thus, $Q_{\mathrm{rot}}$ needs to be divided by $\sigma = 2$. For heteronuclear diatomic molecules you simply have $\sigma = 1$.
With the introduction of $\sigma$ the equation for $Q_{\mathrm{rot}}$ becomes:

\begin{align}
  Q_{\mathrm{rot}} &amp;\approx \frac{1}{\sigma} \left( \frac{T}{\Theta_{\mathrm{rot}}} + \frac{1}{3} + \frac{1}{15} \frac{\Theta_{\mathrm{rot}}}{T} + \frac{4}{315} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} \right) \right) \\
  &amp;= \frac{T}{\sigma \Theta_{\mathrm{rot}}} \left(1  + \frac{1}{3} \frac{\Theta_{\mathrm{rot}}}{T} + \frac{1}{15} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{2} + \frac{4}{315} \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{3} + \mathcal{O}\left( \left(\frac{\Theta_{\mathrm{rot}}}{T}\right)^{4} \right) \right)
\end{align}