Bounds on $S = \frac{1}{1001} + \frac{1}{1002}+ \frac{1}{1003}+\dots+\frac{1}{3001}$

Upper bound $\dfrac{7}{6}$

Your question describes a method of breaking the range up into sub-ranges, and then estimating the sum over a sub-range by noting that each term is $\le$ the first term. But we can do much better if we estimate this sum using a linear approximation instead. Then we only need three sub-ranges to establish the upper bound of $\frac76$.

Consider a 'graph' of the function $$f(k)=\frac{1}{k}$$ for integers $k\ge 1$. If we look at any interval $m\le k\le n$ with $1\le m\le n-2$, we can compare it with the straight-line function joining points $(m,\frac{1}{m})$ and $(n,\frac{1}{n})$. We can define this function explicitly if you want, as: $$g(k)=\frac{1}{n-m}\left(\frac{k-m}{n}+\frac{n-k}{m}\right)$$ Then $f(m)=g(m),f(n)=g(n)$; and for $m<k<n,f(k)<g(k)$ (because $f$ is a convex function).

So $$\sum_{k=m}^nf(k)<\sum_{k=m}^ng(k)$$ But the $g(k)$ are in arithmetic progression, so the RHS is half the sum of the first and the last terms multiplied by the number of terms, i.e. $$\frac{n-m+1}{2}\left(\frac{1}{m}+\frac{1}{n}\right)$$

We divide the range into three equal sub-ranges of size $667$. This gives us: $$\begin{align} \sum_{k=1001}^{3001}\frac{1}{k} & =\sum_{k=1001}^{1667}\frac{1}{k}+\sum_{k=1668}^{2334}\frac{1}{k}+\sum_{k=2335}^{3001}\frac{1}{k}\\ & <\frac{667}{2}\left(\frac{1}{1001}+\frac{1}{1667}\right)+\frac{667}{2}\left(\frac{1}{1668}+\frac{1}{2334}\right)+\frac{667}{2}\left(\frac{1}{2335}+\frac{1}{3001}\right)\\ & <\frac{667}{2}\left(\frac{1}{1000.5}+\frac{1}{1667.5}+\frac{1}{1667.5}+\frac{1}{2334.5}+\frac{1}{2334.5}+\frac{1}{3001.5}\right)\\ & =\frac13+\frac15+\frac15+\frac17+\frac17+\frac19\\ & =\frac{356}{315}\\ & <\frac{7}{6}\end{align}$$

where we have used the fact that $\frac{1}{m}+\frac{1}{n}<\frac{1}{m-\frac12}+\frac{1}{n+\frac12}$ if $m<n$, again by convexity of $f(k)$.

Note that if the last term were $\frac{1}{3000}$, we could do this with just two sub-ranges, and we would get exactly $\frac{7}{6}$ as a bound. So perhaps I have missed an easier way.

Lower bound $\dfrac{29}{27}$

This time we use the convexity of the function to estimate the sum over a sub-range from below. The sum of $\frac{1}{k}$ over the range $m\le k\le n$ is greater than the number of elements multiplied by the middle element if the range contains an odd number of elements. So we get $$\begin{align} \sum_{k=1001}^{3001}\frac{1}{k} & =\sum_{k=1001}^{1667}\frac{1}{k}+\sum_{k=1668}^{2334}\frac{1}{k}+\sum_{k=2335}^{3001}\frac{1}{k}\\ & > 667\left(\frac{1}{1334}+\frac{1}{2001}+\frac{1}{2668}\right)\\ & = \frac12+\frac13+\frac14\\ & = \frac{13}{12}\\ & >\frac{29}{27} \end{align}$$

We can find an approximation of the result.

Consider $$\sum_{i=1}^{2001}\frac1{1000+i}=\sum_{i=1}^{3001}\frac1{i}-\sum_{i=1}^{1000}\frac1{i}$$ and we shall use the fact that $$\sum_{i=1}^{n}\frac1{i}=H_n$$

Now, the asymptotics of harmonic numbers $$H_n=\gamma +\log \left({n}\right)+\frac{1}{2 n}-\frac{1}{12 n^2}+O\left(\frac{1}{n^4}\right)$$ Using it twice and you will have sharper bounds.