The above answers do not mention that we can also find closed form expressions by using harmonic sum formulae.
Introduce the sum
$$S_m(x) = \sum_{n\ge 1} \frac{(nx)^{2m-1}}{e^{2\pi nx}-1}$$
so that we are looking to find $S_m(1).$
The sum term is harmonic and may be evaluated by inverting its Mellin transform.
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 = 1, \quad \mu_k = k \quad \text{and} \quad
g(x) = \frac{x^{2m-1}}{e^{2\pi x}-1}.$$
We need the Mellin transform of $1/(e^{2\pi x}-1)$ which is
$$\int_0^\infty \frac{1}{e^{2\pi x}-1} x^{s-1} dx =
\int_0^\infty \frac{e^{-2\pi x}}{1-e^{-2\pi x}} x^{s-1} dx.$$
Expanding the inner sum we obtain
$$\int_0^\infty \sum_{q\ge 1} e^{-2\pi qx} x^{s-1} dx
= \sum_{q\ge 1} \int_0^\infty e^{-2\pi qx} x^{s-1} dx
= \Gamma(s) \sum_{q\ge 1} \frac{1}{(2\pi q)^s}
= \frac{1}{(2\pi)^s} \Gamma(s) \zeta(s).$$
It follows that the Mellin transform $g^*(s)$ of $g(x)$ is given by
$$\frac{1}{(2\pi)^{s+2m-1}} \Gamma(s+2m-1) \zeta(s+2m-1).$$
Therefore the Mellin transform $Q_m(s)$ of $S_m(x)$ is given by
$$\frac{1}{(2\pi)^{s+2m-1}} \Gamma(s+2m-1) \zeta(s)\zeta(s+2m-1).$$
The Mellin inversion integral here is
$$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q_m(s)/x^s ds$$
which we evaluate by shifting it to the left for an expansion about zero.
We need to examine which poles of the gamma function term are canceled by the two zeta function terms. The poles of the gamma function term are at all integers less than or equal to $-(2m-1).$ The plain zeta term cancels all the even ones and it also cancels the one at $-(2m-1)+1$ from the compound zeta term when $m>1.$ The compound zeta term cancels all the odd ones starting at $-(2m-1)-2.$ This leaves only two poles at $s=1$ and at $s=-(2m-1).$
Therefore we have the following residues that contribute:
$$\mathrm{Res}(Q_m(s)/x^s; s=1) =
\frac{1}{(2\pi)^{2m}} \Gamma(2m) \zeta(2m) \frac{1}{x}
\\= \frac{1}{(2\pi)^{2m}} (2m-1)! \frac{(-1)^{m+1} B_{2m} (2\pi)^{2m}}{2\times (2m)!}
\frac{1}{x} =
\frac{(-1)^{m+1} B_{2m}}{2\times 2m} = \frac{(-1)^{m+1} B_{2m}}{4m} \frac{1}{x}$$
and
$$\mathrm{Res}(Q_m(s)/x^s; s=-(2m-1)) = \zeta(-(2m-1)) \zeta(0) x^{2m-1}
\\= -\frac{1}{2} \times -\frac{B_{2m}}{2m} x^{2m-1}= \frac{B_{2m}}{4m} x^{2m-1}.$$
When $m=1$ there is a pole at $s=0$ which contributes
$$ \mathrm{Res}(Q_m(s)/x^s; s=0) = -\frac{1}{4\pi}.$$
We now specialize in the case when $m$ is odd and at least three. Then we observe a fact which is of critical importance, namely that $Q(s)$ is odd on the line $\Re(s) = -(m-1),$ which we now prove.
Recall the formula for $Q(s),$ which was
$$\frac{1}{(2\pi)^{s+2m-1}} \Gamma(s+2m-1) \zeta(s)\zeta(s+2m-1).$$
Now put $s = -(m-1)+it$ to get
$$\frac{1}{(2\pi)^{m+it}} \Gamma(m+it) \zeta(-m+it+1)\zeta(m+it).$$
Using the functional equation of the Riemann zeta function on the first zeta function term this becomes
$$\frac{2}{(2\pi)^{2m}} \Gamma(m+it) \Gamma(m-it) \zeta(m+it) \zeta(m-it)
\cos\left(\frac{1}{2}(m-it)\pi\right).$$
Now using the series representation of the zeta function in the half plane $\Re(s)>1$ it is easy to see that the two zeta function terms are conjugates of each other and hence their product is a real number. More importantly the value stays the same when $t$ changes sign. Similarly the limit definition of the gamma function shows that the product of the two gamma function terms is a real number, which also stays the same when $t$ changes sign. The factor in front no longer depends on $t.$ That leaves
$$\cos\left(\frac{1}{2}(m-it)\pi\right)
= \cos\left(\frac{\pi}{2}m - \frac{\pi}{2}it\right)
= (-1)^{\frac{m-1}{2}} \sin\left(\frac{\pi}{2}it\right).$$
This term is indeed odd in $t$ which concludes the argument.
Returning to the original computation we now know to shift the integral from $\Re(s) = 3/2$ to $\Re(s) = -(m-1),$ picking up only the residue at $s=1$ and profiting from the fact that the contribution from the left vertical segment is zero. (Set $x=1$ before the shift.) We have established that for $m$ odd and $m>1,$
$$\sum_{n\ge 1} \frac{n^{2m-1}}{e^{2\pi n}-1} = (-1)^{m+1} \frac{B_{2m}}{4m}
= \frac{B_{2m}}{4m}.$$
There is a correction term of $1/\pi/8$ when $m=1$ originating as a half contribution from the pole at zero, which is on the left side of the contour.
$\sum_{0\leq k\leq n/2}\binom{n}{k}x^k = (x+1)^n - \binom{n}{\left\lfloor\frac{n}{2}\right\rfloor + 1} x^{\left\lfloor\frac{n}{2}\right\rfloor + 1} {_2F_1}\left(1,\left\lfloor\frac{n}{2}\right\rfloor - n + 1;\left\lfloor \frac{n}{2}\right\rfloor + 2;-x\right)$,
where $_2F_1$ is the hypergeometric function.
Best Answer
From the series expansion for $\arcsin$, we have
$$\frac{\arcsin(x)}x = \sum_{n=0}^\infty \frac{\binom{2n}n}{2n+1} \left(\frac x2\right)^{2n}$$
and from its square,
$$2\arcsin^2(x) = \sum_{n=0}^\infty \frac{1}{n^2\binom{2n}n} \left(2x\right)^{2n}$$
Both series converge for $x\in[-1,1]$; pick the right value(s) to recover the sums in question.