Q:

How to determine bounds on one variable in a system of inequalities?

I came across this problem whilst exploring the asymptotic behaviour (or not) of different generalised harmonic numbers. I am interested in the point of 'cross-over' between a generalised harmonic number where the denominator of the summand is raised to a power, and a non-exponential harmonic sum operating on some subset of the natural numbers.
For example, take the generalised harmonic number $H_x^{(k)}=\sum_{n=1}^x \frac{1}{n^k}$, and a harmonic number operating only on odd denominators $G_x=\sum_{n=1}^x \frac{1}{2 n-1}$.
Clearly, there exist values of $x,k$ such that $G_x<H_x^{(k)}$ and values such that $H_x^{(k)}<G_x$. Thus there exists a value $c=G_{x_0}$ such that
$$G_{x_0}=c<H_{x_0}^{(k)}=\sum_{n=1}^x \frac{1}{n^k}$$
and
$$H_{{x_0}+2}^{(k)}<G_{{x_0}+2}=c+\frac{1}{2x_0+1}+\frac{1}{2x_0+3}$$
or
$$H_{{x_0}+2}^{(k)}-c<\frac{1}{2x_0+1}+\frac{1}{2x_0+3}$$
The values of $c,x_0,k$ are obviously co-dependent. I am searching for a way to solve for $x_0$ or at least put bounds on it.
I am interested in how to approach this algebraically rather than numerically. This is a single simple example of $G$ and I want to be able to explore how to solve such problems generally, for whatever pattern of $G$ I choose (provided it's formulable!).
Algebraically, how do I put bounds on $x_0$ in terms of $c,k$?
(I am an amateur, so I need a fair amount of hand-holding, hence the bounty.)

A:

I am interested in how to approach this algebraically rather than numerically. This is a single simple example of $G$ and I want to be able to explore how to solve such problems generally, for whatever pattern of $G$ I choose (provided it's formulable!).

We can approximate the sums by the integrals and then to deal with the resulting functions. For instance, for $k<1$,  $$H^{(k)}(x)\approx \int_{1}^{x} \frac 1{t^k} dt=\frac{t^{1-k}}{1-k}{\Huge |}^{x}_1=\frac{1}{1-k}\left(x^{1-k}-1\right).$$
For $k=1$ and natural $x$, according to Wikipedia,
$$H^{(1)}(x)\sim\ln x+\gamma+\frac{1}{2x}- \frac{1}{12x^2}+\frac{1}{120x^4}-\dots,$$
where  $\gamma\approx 0.5772156649$ is the Euler–Mascheroni constant.
For $k>1$ when $x$ tends to infinity, the sequence $H^{(k)}_x$ converges to Riemann zeta function $\zeta(k)$.
Similarly we have $$G(x)\approx \int_1^{x}\frac {1}{2t-1}dt=\frac 12\ln (2x-1).$$

