I'm investigating the behavior of the following function as $x\to \infty$: $$f(x):=\sum_{1\le n\le x}\frac{J_2(n)}{n^2}$$ where $J_k(n)$ is the Jordan's totient function $$J_k(n):=n^k\prod_{p\mid n}(1-\frac1{p^k})$$ and I found that
$$f(x)=\frac x{\zeta (2)}(\ln x+2\gamma-1)+O(\sqrt x)$$
Here's my work:
Since
$$\sum_{d\mid n}J_2(n)=n^2$$
The Möbius inversion formula gives
$$J_2(n)=\sum_{d\mid n}\mu(d) \left(\frac nd\right)^2\ \
\implies\ \ \frac{J_2(n)}{n^2}=\sum_{d\mid n}\frac{\mu(d)}{d^2}$$
Therefore
$$\sum_{n\le x}\frac{J_2(n)}{n^2}=\sum_{n\le x}\sum_{d\mid n}\frac{\mu(d)}{d^2}=\sum_{d\le x}\frac{\mu(d)}{d^2}\sum_{n\le x,d\mid n}1$$
From p44 in the book, the Dirichlet hyperbola method gives
$$\sum_{n\le x,d\mid n}1=x(\ln x+2\gamma -1) +O(\sqrt x)$$
and since
$$\sum_{d\ge 1}\frac{\mu(d)}{d^2}=\frac 1{\zeta(2)}$$
we get the result.
However, when I plot the graph of $$g(x):=\frac{\zeta(2)f(x)-x\ln x}{(2\gamma-1)x}$$ with Mathematica, $g$ doesn't tend to $1$ as $x\to\infty$:
Where did I do wrong?

You should plug in
$$ \sum_{\substack{n\le x\\d|n}}1=\left\lfloor\frac xd\right\rfloor=\frac xd+O(1) $$
instead of the asymptotic formula for the sum over divisor function. This indicates the correct asymptotic expansion should be
$$ \sum_{n\le x}{J_2(n)\over n^2}=x\sum_{d\le x}{\mu(d)\over d^3}+O\left(\sum_{d\ge1}{1\over d^2}\right)={x\over\zeta(3)}+O(1). $$