Sum the series $\sum_{n = 1}^{\infty}\{\coth (n\pi x) + x^{2}\coth(n\pi/x)\}/n^{3}$
It seems to have escaped attention that this sum may be evaluated using harmonic summation techniques.
Put $$S(x) = \zeta(3) + \sum_{n\ge 1} \frac{-1+\coth(n\pi x)}{n^3}$$ and introduce the sum $$T(x) = \sum_{n\ge 1} \frac{-1+\coth(n\pi x)}{n^3}.$$
The sum term is harmonic and may be evaluated by inverting its Mellin transform. We will construct a functional equation for $T(x).$
Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have $$\lambda_k = \frac{1}{k^3}, \quad \mu_k = k \quad \text{and} \quad g(x) = 2\frac{e^{-2\pi x}}{1-e^{-2\pi x}}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$2 \int_0^\infty \frac{e^{-2\pi x}}{1-e^{-2\pi x}} x^{s-1} dx \\ = 2 \int_0^\infty \sum_{q\ge 1} e^{-2q\pi x} x^{s-1} dx = 2 \sum_{q\ge 1} \int_0^\infty e^{-2q\pi x} x^{s-1} dx \\= 2 \Gamma(s) \sum_{q\ge 1} \frac{1}{(2\pi q)^s} = \frac{2}{2^s} \frac{1}{\pi^s} \Gamma(s) \zeta(s).$$
It follows that the Mellin transform $Q(s)$ of the harmonic sum $T(x)$ is given by
$$Q(s) = \frac{2}{2^s} \frac{1}{\pi^s} \Gamma(s) \zeta(s) \zeta(s+3) \quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} \frac{1}{k^3} \frac{1}{k^s} = \zeta(s+3)$$ for $\Re(s) > -2.$
The Mellin inversion integral here is $$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$ which we evaluate by shifting it to the left for an expansion about zero.
Fortunately the trivial zeros of the two zeta function terms cancel the poles of the gamma function term. Shifting to $\Re(s) = -3 -1/2$ we get $$T(x) = \frac{\pi^3 x^3}{90} + 4\zeta'(-2)\pi^2 x^2 + \frac{\pi^3 x}{18} - \zeta(3) + \frac{\pi^3}{90x} + \frac{1}{2\pi i} \int_{-7/2-i\infty}^{-7/2+i\infty} Q(s)/x^s ds.$$
Substitute $s = -2 - t$ in the remainder integral to get $$- \frac{1}{2\pi i} \int_{3/2+i\infty}^{3/2-i\infty} \frac{2}{2^{-2-t}} \frac{1}{\pi^{-2-t}} \Gamma(-2-t) \zeta(-2-t) \zeta(1-t) x^{t+2} dt$$ which is $$\frac{x^2}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} 2^{3+t} \pi^{2+t} \Gamma(-2-t) \zeta(-2-t) \zeta(1-t) x^t dt.$$
In view of the desired functional equation we now use the functional equation of the Riemann zeta function on $Q(s)$ to prove that the integrand of the last integral is in fact $-Q(t).$
Start with the functional equation $$\zeta(1-s) = \frac{2}{2^s\pi^s} \cos\left(\frac{\pi s}{2}\right) \Gamma(s) \zeta(s)$$ and substitute this into $Q(s)$ to obtain $$Q(s) = \frac{2}{2^s} \frac{1}{\pi^s} \frac{\zeta(1-s) 2^s \pi^s}{2\cos\left(\frac{\pi s}{2}\right)} \zeta(s+3) = \frac{\zeta(3+s)}{\cos\left(\frac{\pi s}{2}\right)} \zeta(1-s).$$ Apply the functional equation again (this time to $\zeta(s+3)$) to get $$Q(s) = \frac{1}{\cos\left(\frac{\pi s}{2}\right)} \frac{2}{2^{-2-s} \pi^{-2-s}} \cos\left(\frac{\pi (-2-s)}{2}\right) \Gamma(-2-s) \zeta(-2-s) \zeta(1-s)$$ Observe that $$\frac{\cos\left(-\pi-\frac{\pi s}{2}\right)} {\cos\left(\frac{\pi s}{2}\right)} = - \frac{\cos\left(-\frac{\pi s}{2}\right)} {\cos\left(\frac{\pi s}{2}\right)} = -1$$ so we finally get $$Q(s) = - 2^{3+s} \pi^{2+s} \Gamma(-2-s) \zeta(-2-s) \zeta(1-s),$$ thus proving the claim.
We have established the functional equation $$T(x) = \frac{\pi^3 x^3}{90} + 4\zeta'(-2)\pi^2 x^2 + \frac{\pi^3 x}{18} - \zeta(3) + \frac{\pi^3}{90x} - x^2 T(1/x).$$
Finally returning to the sum that was the initial goal we see that it has the value $$\zeta(3) + T(x) + x^2 (\zeta(3) + T(1/x))$$ or $$\zeta(3) + T(x) + x^2 \zeta(3) + x^2 T(1/x).$$ Using the functional equation for $T(x)$ this becomes $$\zeta(3) + T(x) + x^2 \zeta(3) + \frac{\pi^3 x^3}{90} + 4\zeta'(-2)\pi^2 x^2 + \frac{\pi^3 x}{18} - \zeta(3) + \frac{\pi^3}{90x} - T(x)$$ which is $$x^2 \zeta(3) + \frac{\pi^3 x^3}{90} + 4\zeta'(-2)\pi^2 x^2 + \frac{\pi^3 x}{18} + \frac{\pi^3}{90x}.$$
The inspiration for this calculation is from the paper "Mellin Transform and its Applications" by Szpankowski.
Addendum. In view of the fact that $$\zeta(3) + 4\zeta'(-2)\pi^2 =0 $$ (consult e.g. MathWorld) this finally becomes $$\frac{\pi^3 x^3}{90} + \frac{\pi^3 x}{18} + \frac{\pi^3}{90x} = \frac{\pi^3}{90x} \left(x^4 + 5x^2 + 1\right).$$
Addendum II. There is another functional equation of a harmonic sum at this MSE link, this one somewhat more advanced.
Following in the same manner as this answer...
We are going to use the contour integral $$ \oint\pi\cot\left(\frac{\pi z}{\pi x}\right)\left(\frac{\coth(z)}{z^3}-\frac1{z^4}-\frac1{3z^2}\right)\mathrm{d}z=0\tag{1} $$ where the contours of interest are, for real $R\to\infty$ and integer $n\to\infty$, $$ \small\textstyle\color{#00A000}{[R,-R]+(n+\frac12)\pi i}\cup\color{#C00000}{-R+(n+\frac12)\pi i[1,-1]}\cup\color{#00A000}{[-R,R]-(n+\frac12)\pi i}\cup\color{#C00000}{R+(n+\frac12)\pi i[-1,1]} $$ The integral along the red paths becomes negligible as $R\to\infty$. Along the upper green path, where $\mathrm{Im}(z)\approx+\infty$, $\cot(z)\approx-i$. Along the lower green path, where $\mathrm{Im}(z)\approx-\infty$, $\cot(z)\approx+i$. Since $\coth(z+\frac\pi2i)=\tanh(z)$, the integral along each of the green paths tends to $0$. Therefore, the full integral is $0$.
Since $$ \pi\cot\left(\frac{\pi z}{\pi x}\right)\text{ has residue }\pi x\text{ at }z=\pi nx\tag{2} $$ and $$ \frac{\coth(z)}{z^3}-\frac1{z^4}-\frac1{3z^2}=-\frac1{45}+O(z^2)\text{ at }z=0\tag{3} $$ the contribution from the singularities on the real axis is $$ 2\pi i\cdot\pi x\left[2\sum_{n=1}^\infty\left(\frac{\coth(\pi nx)}{(\pi nx)^3}-\frac1{(\pi nx)^4}-\frac1{3(\pi nx)^2}\right)-\frac1{45}\right]\tag{4} $$ Since $$ \frac{\coth(z)}{z^3}\text{ has residue }\frac1{(\pi in)^3}\text{ at }z=\pi i n\text{ for }n\ne0\tag{5} $$ and $$ \pi\cot\left(\frac{\pi z}{\pi x}\right)=-\pi i\coth\left(\frac{\pi n}{x}\right)\text{ at }z=\pi in\tag{6} $$ the contribution from the singularities on the imaginary axis is $$ 2\pi i\left[2\sum_{n=1}^\infty\frac\pi{x^3}\frac{\coth\left(\frac{\pi n}{x}\right)}{\left(\frac{\pi n}{x}\right)^3}\right]\tag{7} $$ Combining $(1)$, $(4)$, and $(7)$, yields $$ x^2\sum_{n=1}^\infty\frac{\coth(\pi nx)}{(\pi nx)^3}+\frac1{x^2}\sum_{n=1}^\infty\frac{\coth\left(\frac{\pi n}{x}\right)}{\left(\frac{\pi n}{x}\right)^3} =\frac{\zeta(4)}{\pi^4x^2}+\frac{\zeta(2)}{3\pi^2}+\frac{x^2}{90}\tag{8} $$ Multiplying by $\pi^3x$ to match the question, we get $$ \sum_{n=1}^\infty\frac{\coth(\pi nx)+x^2\coth(\pi n/x)}{n^3}=\frac{\pi^3}{90x}\left(1+5x^2+x^4\right)\tag{9} $$
Recall the well known Mittag-Leffler expansion of hyperbolic cotangent function (denote $\mathbb{W}=\mathbb{Z}/\{0\}$) :
$$\sum_{m\in\mathbb{W}}\frac{1}{m^2+z^2}=\frac{\pi\coth\pi z}{z}-\frac{1}{z^2}\tag{ML}$$
Hence, your sum is by its symmetry :
$$\begin{align} S&=\frac{1}{2}\sum_{n \in \mathbb{W}}\{\coth (n\pi x) + x^{2}\coth(n\pi/x)\}/n^{3} \\ \\ &=\frac{1}{2\pi x}\sum_{n \in \mathbb{W}}\left(\frac{1}{n^4}+\sum_{m\in\mathbb{W}}\frac{x^2/n^2}{m^2+n^2x^2}\right) +\left(\frac{x^4}{n^4}+\sum_{m\in\mathbb{W}}\frac{x^2/n^2}{m^2+n^2/x^2}\right)\tag{1}\\ \\ &=\frac{1}{2\pi x}\left(\zeta(4)+x^4\zeta(4)+\sum_{n,m \in \mathbb{W}^2}\frac{x^2}{n^2}\frac{1}{m^2+n^2x^2}+\frac{x^2}{n^2}\frac{1}{m^2+n^2/x^2}\right)\tag{2}\\ \\ &=\frac{1}{2\pi x}\left(2\zeta(4)+2x^4\zeta(4)+x^2\sum_{n,m\in\mathbb{W}^2}\frac{1}{n^2m^2}\frac{m^2+n^2x^2}{m^2+n^2x^2}\right)\tag{3}\\ \\ &=\frac{1}{2\pi x}\left(2\zeta(4)+2x^4\zeta(4)+4x^2\zeta^2(2)\right)\tag{4}\\ \\ &=\frac{1}{2\pi x}\left(2\frac{\pi^4}{90}+2x^4\frac{\pi^4}{90}+4x^2\frac{\pi^4}{36}\right)\tag{5}\\ \\ &=\frac{\pi^3}{90x}\left(1+x^4+5x^2\right) \end{align}$$
Explanations
$(1)$ Use the Mittag-Leffler formula (ML) with $z=nx$ and $z=n/x$
$(2,4)$ Recall $\zeta(s)=\sum_{n=1}^{\infty}1/n^s$
$(3)$ In the second sum rename $n \longleftrightarrow m$
$(5)$ Zetas for $s=2$ and $4$ are well known, i.e. $\zeta(2)=\pi^2/6$ and $\zeta(4)=\pi^4/90$