Ramanujan found the following formula:
$$\large \sum_{n=1}^\infty \frac{n^{13}}{e^{2\pi n}-1}=\frac 1{24}$$
I let $e^{2\pi n}-1=\left(e^{\pi n}+1\right)\left(e^{\pi n}-1\right)$ to try partial fraction decomposition and turn the sum into telescoping, but methinks it doesn't lead anywhere and only makes things hairy.
How does one go about proving this? Thanks.
Suppose we seek to evaluate
$$S = \sum_{n\ge 1} \frac{n^{13}}{e^{2\pi n}-1}.$$
This sum may be evaluated using harmonic summation techniques.
Introduce the sum $$S(x; p) = \sum_{n\ge 1} \frac{n^{4p+1}}{e^{nx}-1}$$ with $p$ a positive integer and $x\gt 0.$
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 = k^{4p+1}, \quad \mu_k = k \quad \text{and} \quad g(x) = \frac{1}{e^x-1}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty \frac{1}{e^{x}-1} x^{s-1} dx = \int_0^\infty \frac{e^{-x}}{1-e^{-x}} x^{s-1} dx \\ = \int_0^\infty \sum_{q\ge 1} e^{-q x} x^{s-1} dx = \sum_{q\ge 1} \int_0^\infty e^{-q x} x^{s-1} dx \\= \Gamma(s) \sum_{q\ge 1} \frac{1}{q^s} = \Gamma(s) \zeta(s).$$
It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x,p)$ is given by
$$Q(s) = \Gamma(s) \zeta(s) \zeta(s-(4p+1)) \\ \text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} k^{4p+1} \frac{1}{k^s} = \zeta(s-(4p+1))$$ for $\Re(s) > 4p+2.$
The Mellin inversion integral here is $$\frac{1}{2\pi i} \int_{4p+5/2-i\infty}^{4p+5/2+i\infty} Q(s)/x^s ds$$ which we evaluate by shifting it to the left for an expansion about zero.
The two zeta function terms cancel the poles of the gamma function term and we are left with just
$$\begin{align} \mathrm{Res}(Q(s)/x^s; s=4p+2) & = \Gamma(4p+2) \zeta(4p+2) / x^{4p+2} \quad\text{and}\\ \mathrm{Res}(Q(s)/x^s; s=0) & = \zeta(0) \zeta(-(4p+1)). \end{align}$$
Computing these residues we get
$$(4p+1)! \frac{B_{4p+2} (2\pi)^{4p+2}}{2(4p+2)! \times x^{4p+2}} = \frac{B_{4p+2} (2\pi)^{4p+2}}{2\times (4p+2) \times x^{4p+2}}$$ and $$- \frac{1}{2} \times -\frac{B_{4p+2}}{4p+2}.$$
This shows that $$S(x;p) = \frac{B_{4p+2} (2\pi)^{4p+2}}{(8p+4)\times x^{4p+2}} + \frac{B_{4p+2}}{8p+4} + \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} Q(s)/x^s ds.$$
To treat the integral recall the duplication formula of the gamma function: $$\Gamma(s) = \frac{1}{\sqrt\pi} 2^{s-1} \Gamma\left(\frac{s}{2}\right) \Gamma\left(\frac{s+1}{2}\right).$$
which yields for $Q(s)$
$$\frac{1}{\sqrt\pi} 2^{s-1} \Gamma\left(\frac{s}{2}\right) \Gamma\left(\frac{s+1}{2}\right) \zeta(s) \zeta(s-(4p+1))$$
Furthermore observe the following variant of the functional equation of the Riemann zeta function: $$\Gamma\left(\frac{s}{2}\right)\zeta(s) = \pi^{s-1/2} \Gamma\left(\frac{1-s}{2}\right) \zeta(1-s)$$
which gives for $Q(s)$ $$\frac{1}{\sqrt\pi} 2^{s-1} \pi^{s-1/2} \Gamma\left(\frac{s+1}{2}\right) \Gamma\left(\frac{1-s}{2}\right) \zeta(1-s)\zeta(s-(4p+1)) \\ = \frac{1}{\sqrt\pi} 2^{s-1} \pi^{s-1/2} \frac{\pi}{\sin(\pi(s+1)/2)} \zeta(1-s)\zeta(s-(4p+1)) \\ = 2^{s-1} \frac{\pi^s}{\sin(\pi(s+1)/2)} \zeta(1-s)\zeta(s-(4p+1)).$$
Now put $s=4p+2-u$ in the remainder integral to get
$$- \frac{1}{x^{4p+2}} \frac{1}{2\pi i} \int_{4p+5/2+i\infty}^{4p+5/2-i\infty} 2^{4p+1-u} \\ \times \frac{\pi^{4p+2-u}}{\sin(\pi(4p+3-u)/2)} \zeta(u-(4p+1))\zeta(1-u) x^u du \\ = \frac{2^{4p+2} \pi^{4p+2}}{x^{4p+2}} \frac{1}{2\pi i} \int_{4p+5/2-i\infty}^{4p+5/2+i\infty} 2^{u-1} \\ \times \frac{\pi^{u}}{\sin(\pi(4p+3-u)/2)} \zeta(u-(4p+1))\zeta(1-u) (x/\pi^2/2^2)^u du.$$
Now $$\sin(\pi(4p+3-u)/2) = \sin(\pi(1-u)/2+\pi (2p+1)) \\ = - \sin(\pi(1-u)/2) = \sin(\pi(-1-u)/2) = - \sin(\pi(u+1)/2).$$
We have shown that $$\bbox[5px,border:2px solid #00A000] {S(x;p) = \frac{B_{4p+2} (2\pi)^{4p+2}}{(8p+4)\times x^{4p+2}} + \frac{B_{4p+2}}{8p+4} - \frac{(2\pi)^{4p+2}}{x^{4p+2}} S(4\pi^2/x;p)}.$$
In particular we get
$$S(2\pi; p) = \frac{B_{4p+2}}{8p+4}.$$
The sequence in $p$ starting from $p=1$ is
$${\frac{1}{504}},{\frac{1}{264}},1/24, {\frac{43867}{28728}},{\frac{77683}{552}}, {\frac{657931}{24}},{\frac{1723168255201}{171864}}, \ldots$$
We thus have for $p=3$ as per request from OP
$$\bbox[5px,border:2px solid #00A000]{ \sum_{n\ge 1} \frac{n^{13}}{e^{2\pi n}-1} = \frac{1}{24}.}$$
References, as per request, are: Flajolet and Sedgewick, Mellin transform asymptotics, INRIA RR 2956 and Szpankowski, Mellin Transform and its applications, from Average Case Analysis of Algorithms on Sequences.
It is the weight $14$ Eisenstein series $$G_{14}(z)=\sum_{(n,m)\ne (0,0)} \frac1{(zn+m)^{14}}= 2\zeta(14)+\sum_{n\ne 0} \frac{1}{13!} \frac{d^{13}}{dz^{13}}\frac{2i\pi}{e^{2i\pi n z}-1}$$ $$=2\zeta(14)+\sum_{n\ge 1} \frac{4i\pi}{13!} \sum_{m\ge 1} (2i\pi m)^{13}e^{2i\pi mz}=2\zeta(14)+(2i\pi)^{14}\frac{2}{13!}\sum_{k\ne 1}\frac{k^{13}}{e^{-2i\pi kz}-1} $$
$$G_{14}(z)= z^{-14}G_{14}(-1/z)\implies \qquad G_{14}(i)=0$$
$$\boxed{(2i\pi)^{14}\frac{2}{13!}\sum_{k\ne 1}\frac{k^{13}}{e^{2\pi kz}-1}+2\zeta(14)=0 }$$ $2\zeta(14)=-\frac{B_{14}(2\pi)^{14}}{(14)!} $
For your curiosity !
I do not know if these results are known but, beside this one, $$ \sum_{n=1}^\infty \frac{n^{5}}{e^{2\pi n}-1}=\frac 1{504}=\frac 1{21 \times 24}\qquad\text{and} \qquad \sum_{n=1}^\infty \frac{n^{9}}{e^{2\pi n}-1}=\frac 1{264}=\frac 1{11 \times 24}$$
If they are known, please tell me where I could find them.
Theorem 1. (see [1] pg.275-276)
Let $a,b>0$ with $ab=\pi^2$, and let $\nu$ be any non zero integer. Then
$$
a^{-\nu}\left\{\frac{1}{2}\zeta(2\nu+1)+\sum^{\infty}_{n=1}\frac{n^{-2\nu-1}}{e^{2an}-1}\right\}-
(-b)^{-\nu}\left\{\frac{1}{2}\zeta(2\nu+1)+\sum^{\infty}_{n=1}\frac{n^{-2\nu-1}}{e^{2bn}-1}\right\}=
$$
\begin{equation}
=-2^{2\nu}\sum^{\nu+1}_{n=0}(-1)^n\frac{B_{2n}}{(2n)!}\frac{B_{2\nu+2-2n}}{(2\nu+2-2n)!}a^{\nu+1-n}b^n,\tag 1
\end{equation}
where $\zeta(s)$ is the Riemann zeta function and $B_n$ is the $n-$th Bernoulli number.
Notes
For integer $\nu<-1$ formula (1) evaluated from:
Theorem 2. (see [1] pg.261) If $\nu$ is integer greater than 1, then ($ab=\pi^2$, $a,b>0$) $$ a^{\nu}\sum^{\infty}_{n=1}\frac{n^{2\nu-1}}{e^{2an}-1}-(-b)^{\nu}\sum^{\infty}_{n=1}\frac{n^{2\nu-1}}{e^{2bn}-1}=(a^{\nu}-(-b)^{\nu})\frac{B_{2\nu}}{4\nu}\tag 2 $$
[1]: B.C. Berndt, 'Ramanujan`s Notebooks Part II'. Springer Verlang, New York., (1989).
Essentially the same as @reuns answer but with more detail.
Let $k\ge2$ be an integer and define the Eisenstein series $$G_{2k}(z)=\sum_{(n,m)\in A}\frac{1}{(n+mz)^{2k}},\tag 1$$ where $A=\Bbb Z^2\setminus\{(0,0)\}$, and $z\in\Bbb C$ with $\Im(z)>0$. It is simple to show that $G_{2k}(z+1)=G_{2k}(z)$ for all $z$, so it follows that we may write $$G_{2k}(z)=\sum_{n\ge0}g_nq^n,$$ where $q=e^{2i\pi z}$. It can be shown (see here) that there is a closed form for $g_n$. Namely, we can write $$\begin{align} G_{2k}(z)&=2\zeta(2k)\left(1+c_{2k}\sum_{n\ge1}\sigma_{2k-1}(n)q^n\right)\\ &=2\zeta(2k)\left(1+c_{2k}\sum_{n\ge1}\frac{n^{2k-1}q^n}{1-q^n}\right), \end{align}$$ where $c_{2k}=\frac{(2\pi i)^{2k}}{(2k-1)!\zeta(2k)}=\frac{-4k}{B_{2k}}=\frac{2}{\zeta(1-2k)}$.
On the other hand, it is simple to show from $(1)$ that $$G_{2k}(-1/z)=z^{2k}G_{2k}(z).$$ Letting $E_{2k}(z)=\frac{1}{2\zeta(2k)}G_{2k}(z)$ for convenience, we have $$E_{2k}(-1/z)=z^{2k}E_{2k}(z).\tag2$$ Then defining $$S_k(q)=\sum_{n\ge1}\frac{n^{2k-1}}{q^n-1},$$ we have $$E_{2k}(z)=1+c_{2k}S_k(e^{-2i\pi z}).$$ Then from $(2)$, we have $$1+c_{2k}S_k(e^{2i\pi/z})=z^{2k}\left(1+c_{2k}S_k(e^{-2i\pi z})\right).\tag3$$ Since your sum is given by $S_7(e^{2\pi})$, we set $k=7$ and $z=i$ in $(3)$, and get $$\begin{align} 1+c_{14}S_7(e^{2\pi})&=-\left(1+c_{14}S_7(e^{2\pi})\right)\\ \Rightarrow S_7(e^{2\pi})&=-\frac{1}{c_{14}}. \end{align}$$ Since $c_{2k}=-4k/B_{2k}$, we have $-1/c_{14}=B_{14}/28=1/24$, and thus $$S_7(e^{2\pi})=\sum_{n\ge1}\frac{n^{13}}{e^{2\pi n}-1}=\frac{1}{24}.$$ In general, the same method allows us to compute $$S_{2k+1}(e^{2\pi})=\sum_{n\ge1}\frac{n^{4k+1}}{e^{2\pi n}-1}=\frac{B_{4k+2}}{8k+4},$$ as was shown by @MarcoRiedel.
This is not a strict solution but a heuristic proof using CAS. It also shows that there are other "magic" numbers instead of 13 which result in similar simple fractions. Playing around I also found the interesting conincidence that in some cases the sum and the related integral give the same results.
Defining
$$S(m) = \sum_{n=1}^{\infty} \frac{n^m}{e^{2 \pi n}-1}\tag{1}$$
we have to show that $S(13) = \frac{1}{24}$.
My first idea was to identify the denominator. I knew this expression from the Planck radioation formula and from the generating function of the Benoulli numbers.
But why so complicated? It is just the sum of a power series. Indeed we can write
$$\frac{1}{e^{2 \pi n}-1} = \sum _{j=1}^{\infty } \exp (- 2 \pi j n)\tag{2}$$
Next, for the numerator we replace the power of $n$ by $z^n$ and consider the intermediate generating function
$$g_{0}(z,j) = \sum _{n=1}^{\infty } z^n \exp (- 2 \pi j n)=\frac{z}{e^{2 \pi j}-z}\tag{3}$$
Summing over $j$ gives the generating function
$$\begin{align}g(z) = & \sum _{j=1}^{\infty } \frac{z}{e^{2 \pi j}-z}\\=& -\frac{1}{2 \pi}\psi _{e^{-2 \pi }}^{(0)}\left(1-\frac{\log (z)}{2 \pi }\right)+\frac{1}{2 \pi} \log \left(1-e^{-2 \pi }\right) \end{align}\tag{4}$$
Here $\psi _{q }^{(0)}(x)$ is the q-digamma function.
The the sums in question here can be found as derivatives of $g(z)$
$$S(m,z) = (z\frac{\partial}{\partial z})^m g(z) |_{z \to 1}\tag{5}$$
Using Mathematica for some values of $m$ I saw the pattern and found the general formula
$$S(m,z)=\frac{(-1)^{m+1} \psi _{e^{-2 \pi }}^{(m)}\left(1-\frac{\log (z)}{2 \pi }\right)}{2^{m+1} \pi ^{m+1}}\tag{6}$$
which gives at $z=1$ where the $\log$ vanishes the result
$$S(m)=S(m,z\to 1) = \frac{1}{(2 \pi )^{n+1}} \psi _{e^{-2 \pi }}^{(n)}(1)\tag{7}$$
Finally, the numerical results in Mathematica for $\frac{1}{S(m)}$ were particularly simple for 3 values of $m$. In the format {m,1/S(m)} we have
$$\{\{5, 504\}, \{9, 264\}, \{13,24\}\}\tag{8}$$
The case $m=13$ leads to the desired result. The "magic" values of $m$ have the form $1+4k$.
Discussion
It is always tempting with sums to look at the corresponding integral.
In our case we consider
$$i(m) = \int_0^{\infty } \frac{n^m}{\exp (2 \pi n)-1} \, dn\tag{9}$$
where in contrast to the sum the integration starts at $n=0$.
The surprising observation is that we have (checked numerically)
$$i(m) = S(m),m=1+4k, k=1,2,...\tag{10}$$
which includes the "magic" values.
It would be nice to find a proof of this observation
