<s>[INST] Deriving a Series Representation of an Integral (Riemann Zeta Function) It is known that
$$\int_0^\infty\frac{x^n}{e^x-1}dx=n!\zeta(n+1)$$
for integer values of $n$ (this is also generalises, but that's not important for this question). I have also had a look at how to arrive at this expression, starting from the series representation of the Riemann Zeta function.
However, just out of interest, I've tried to see if I could manage to derive the series representation for the case that sparked my interest in the first place $(n=3)$ on my own - that is, working backwards from this formula to the series that initially defines the Riemann Zeta function. I did manage to arrive at a series representation, but it's not the correct series.
Here's what I did:
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\int_0^\infty\frac{x^3(e^x+1)}{e^{2x}-1}dx\\
&=\int_0^\infty\frac{x^3}{e^{2x}-1}dx+\int_0^\infty\frac{x^3e^x}{e^{2x}-1}dx\\
&=\frac{1}{16}\int_0^\infty\frac{x^3}{e^x-1}dx+\int_0^\infty\frac{x^3e^x}{e^{2x}-1}dx
\end{align}
Subtracting the first term on the RHS, we have
\begin{align}
\frac{15}{16}\int_0^\infty\frac{x^3}{e^x-1}dx=\int_0^\infty\frac{x^3e^x}{e^{2x}-1}dx
\end{align}
Multiplying across the fraction, we have
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\frac{16}{15}\int_0^\infty\frac{x^3e^x}{e^{2x}-1}dx\\
&=\frac{1}{15}\int_0^\infty\frac{x^3e^{\frac{x}{2}}}{e^x-1}dx
\end{align}
Letting $z:=e^x$, then $x=\log z\implies dx=\frac{dz}{z}$.
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\frac{1}{15}\int_1^\infty\frac{\log^3z\sqrt{z}}{z(z-1)}dz
\end{align}
Letting $u:=\frac{1}{z}$, we have $z=\frac{1}{u}\implies dz=-\frac{du}{u^2}$, $u(1)=1$ and $u(z\rightarrow\infty)=0$. Thus
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\frac{1}{15}\int_0^1\frac{\log^3\left(\frac{1}{u}\right)\sqrt{u}}{u^2}\frac{1}{1-\frac{1}{u}}du\\
&=-\frac{1}{15}\int_0^1\frac{\log^3u\sqrt{u}}{u}\frac{1}{u-1}du\\
&=\frac{1}{15}\int_0^1\frac{\log^3u}{\sqrt{u}}\frac{1}{1-u}du\\
&=\frac{1}{15}\int_0^1\frac{\log^3u}{\sqrt{u}}\sum_{k=0}^\infty u^k du\\
&=\frac{1}{15}\sum_{k=0}^\infty\int_0^1u^{k-\frac{1}{2}}\log^3u du
\end{align}
Letting $t:=\log u$, then $u=e^t\implies du=e^tdt$, $t(u\rightarrow 0)=-\infty$ and $t(1)=0$.
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\frac{1}{15}\sum_{k=0}^\infty\int_{-\infty}^0 t^3e^{\left(k+\frac{1}{2}\right)t}dt
\end{align}
Letting $p:=-t$, then $dt=-dp$, $p(t\rightarrow -\infty)=\infty$ and $p(0)=0$.
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\frac{1}{15}\sum_{k=0}^\infty\int_0^\infty p^3e^{-\left(k+\frac{1}{2}\right)p}dp
\end{align}
Letting $s:=\left(k+\frac{1}{2}\right)p$, then $dp=\frac{ds}{k+\frac{1}{2}}$. The bounds of integration are unchanged, and thus
\begin{align}
\int_0^\infty\frac{x^3}{e^x-1}dx&=\frac{1}{15}\sum_{k=0}^\infty\frac{1}{\left(k+\frac{1}{2}\right)^4}\int_0^\infty s^3e^{-s}ds\\
&=\frac{\Gamma(4)}{15}\sum_{k=0}^\infty\frac{1}{\left(k+\frac{1}{2}\right)^4}\\
&=\frac{2}{5}\sum_{k=0}^\infty\frac{1}{\left(k+\frac{1}{2}\right)^4}
\end{align}
EDIT: Apparently, this is simply a convoluted way of deriving the correct answer, as the series evaluates to $\frac{\pi^4}{6}$, making the final expression equal to the known value of
$$\int_0^\infty\frac{x^3}{e^x-1}dx = 6\sum_{k=1}^\infty\frac{1}{k^4} = \frac{\pi^4}{15}$$
Many thanks to Simple Art for verifying this result!
 [/INST] Your so called proof is very much unorthodox.  A much shorter proof can be found in this post and it basically goes as follows:
$$\begin{align}\int_0^\infty\frac{x^k}{e^x-1}\ dx&=\int_0^\infty\frac{x^ke^{-x}}{1-e^{-x}}\ dx\\&=\int_0^\infty x^ke^{-x}\sum_{n=0}^\infty e^{-nx}\ dx\\&=\int_0^\infty\sum_{n=1}^\infty x^ke^{-nx}\ dx\\&=\int_0^\infty\sum_{n=1}^\infty\frac1{n^{s+1}}y^{k+1}e^{-y}\ dy\tag{$y=nx$}\\&=\sum_{n=1}^\infty\frac1{n^{s+1}}\int_0^\infty y^{k+1}e^{-y}\ dy\\&=\zeta(s+1)\Gamma(s+1)\end{align}$$
The last line comes by definition.  And when $s\in\mathbb N$, $\Gamma(s+1)=s!$
Around the middle part, where you do substitutions, there appears to be missing factors (see comments).
In essence, I think the above is a much simpler derivation that uses the exact same geometric expansion, just in a far simpler way.

As you have edited your question, I give you the relief that your solution is correct.  Particularly, scroll down quite a ways and you will find a nice surprise.
</s>