I've copy the identity in my Question from the solution of Wolfram Alpha online calculator. The expression is tedious to write thus I hope that there are no typos. When you type the code
sum 1/((k+x)(k^2+x)) for k=1 to infinite
and after you can copy the solution from the online calculator with the option $\text{Plaintext}$, you can compute the limit as
lim x--->0 (-sqrt(-x) polygamma(0, 1-sqrt(-x))-polygamma(0, 1-sqrt(-x))+sqrt(-x) polygamma(0, sqrt(-x)+1)-polygamma(0, sqrt(-x)+1)+2 polygamma(0, x+1))/(2 (sqrt(-x)-1) (sqrt(-x)+1) x)
Thus, from LHS, one deduces $$\lim_{x\to 0}\sum_{k=1}^\infty\frac{1}{(k+x)(k^2+x)}=\zeta(3).$$
Question. Let $\psi^{(n)}(x)$ the nth derivative of the digamma function. Please can you explain how one deduces rigurously $$\lim_{x\to 0}\frac{-(1+\sqrt{-x})\psi^{(0)}(1-\sqrt{-x})+(1+\sqrt{-x})\psi^{(0)}(\sqrt{-x}+1)+2\psi^{(0)}(1+x)}{2(\sqrt{-x}-1)(\sqrt{-x}+1)x}=-\frac{\psi^{(2)}(1)}{2}$$ holds? Thanks in advance.
Thus in previous question I am asking how you solves $\frac{0}{0}$. I presume that is using L'Hôpital's rule and your knowledges of this specials functions. And I know that $\psi^{(2)}(1)=-2\zeta(3)$, thus I am asking about how compute previous limit rigurously: please can you explain us why are continuous those functions in the limit, each summand, product and quotient, y what is the limit that you are computing (I know that $\lim_{x\to 0^+}=\lim_{x\to 0^-}$ from LHS since LHS is continuous) that you can show to provide us a rigurously proof of the identity.
Just for your curiosity.
As you wrote, if $$S=\sum_{k=1}^\infty\frac{1}{(k+x)(k^2+x)}$$ the expression you give can simplify as $$S=\frac{\left(1+\sqrt{-x}\right) H_{-\sqrt{-x}}+\left(1-\sqrt{-x}\right) H_{\sqrt{-x}}-2 H_x}{2 x (x+1)}$$ using generalized harmonic numbers.
Now, using Taylor expansion around $x=0$ (a CAS did it) you have $$S=-\frac{\psi ^{(2)}(1)}{2}+x \left(\frac{\psi ^{(4)}(1)}{24}-\frac{\pi ^4}{90}\right)+O\left(x^{3/2}\right)$$ which can also write $$S=\zeta (3)- \left(\zeta (4)+\zeta (5)\right)x+O\left(x^{3/2}\right)$$ What is amazing is that, increasing the order of the expansion, we should get $$S=\zeta (3)- \left(\zeta (4)+\zeta (5)\right)x+\left(\zeta (5)+\zeta (6)+\zeta (7)\right)x^2 +O\left(x^{5/2}\right)$$
Edit
This can be generalized to $$S=\sum_{k=1}^\infty\frac{1}{(k^a+x)(k^b+x)}$$ $$S=\zeta (a+b)- \left[\zeta (2 a+b)+\zeta (a+2 b)\right]x+ \left[\zeta (2a+2b)+\zeta (3 a+b)+\zeta (a+3 b)\right]x^2+O\left(x^{3}\right)$$ which is easy to obtain if we start developing $\frac{1}{(k^a+x)(k^b+x)}$ as a Taylor series around $x=0$.