The integral : $\frac{1}{2}\int_0^\infty x^n \operatorname{sech}(x)\mathrm dx$ How can I evaluate
$$\frac{1}{2}\int_0^\infty x^n \operatorname{sech}(x)\mathrm dx?$$
I was trying integration by parts but it seemed like it is getting more complicated.
$$\int_0^\infty x^n \operatorname{sech}(x)\mathrm dx=\left.2x^n\arctan\left(\tanh(x/2)\right)\right|_0^\infty-2n\int_0^\infty x^{n-1}\arctan(\tanh(x/2))\mathrm dx$$
Herein, it seems like we have to apply integration by parts $n$ times but it is not practically possible.

This question is a more general problem of the integral $\int_0^\infty \frac{x}{e^x+e^{-x}}\mathrm dx$, which I was first solving.
Let me know if there's any other method for evaluating this integral. It will be highly appreciated.

I have posted my solution employing a method using Geometric series to which this Wikipedia article helped me in finding the solution.
Please see my answer below.
 A: Alternatively
\begin{align}
\frac{1}{2}\int_0^\infty \frac{x^n}{\cosh x}dx \overset{t=e^{-x}}=&\int_0^1 \frac{(-1)^n\ln^n t}{1+t^2}dt=n!\>\text{Im}\>\text{Li}_{n+1}(i)
\end{align}
A: \begin{align}\mathcal{I}&=\frac{1}{2}\int_0^\infty \frac{x^n}{\cosh(x)}\mathrm dx\\&=\int_0^\infty \frac{x^n}{e^x+e^{-x}}\mathrm dx\\&=\int_0^\infty \frac{x^ne^{-x}}{1+e^{-2x}}\mathrm dx.\end{align}

Proposition: $$\int_0^\infty \frac{x^{s-1}e^{-x}}{1+e^{-2x}}\mathrm dx=\beta(s)\Gamma(s),$$ where $\beta(s)$ is the Dirichlet beta function defined as $\beta(s)=\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^s}$.
Proof:\begin{align}\int_0^\infty x^{s-1}\left(\frac{e^{-x}}{1-(-e^{-2x})}\right)\mathrm dx&=\int_0^\infty x^{s-1}\left(\sum_{k=0}^\infty (-1)^ke^{-(2k+1)x}\right)\mathrm dx \text{, using geometric series}\\&=\sum_{k=0}^\infty (-1)^k\int_0^\infty x^{s-1} e^{-(2k+1)x}\mathrm dx\\&=\sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)^s}\displaystyle\int_0^\infty x^{s-1}e^{-x}\mathrm dx\text{, substituting $(2k+1)x\mapsto x$}\\&=\beta(s)\Gamma(s). \end{align}

Therefore,
$$\mathcal{I}=\beta(n+1)\Gamma(n+1).$$
We can evaluate for different values of $n\ge 0$.
For $n=1$, $\frac{1}{2}\int_0^\infty x\operatorname{sech}(x)=\beta(2)\Gamma(2)=G$, where $G$ is Catalan's constant.
For $n=2$, $\int_0^\infty x^2\operatorname{sech}(x)=2\beta(3)\Gamma(3)=\frac{\pi^3}{8}$
A: Letting $\tan \theta=e^{-x}$ transforms the integral into
$$
\int_0^{\infty} \frac{x^n}{e^x+e^{-x}}dx=(-1)^n \int_0^{\frac{\pi}{4}} \ln^n (\tan \theta) d \theta
$$
By my post,
$$\int_{0}^{\frac{\pi}{4}} \ln ^{n}(\tan \theta) d \theta =(-1)^{n} \beta(n+1) \Gamma(n+1)$$
Hence$$\boxed{\int_0^{\infty} \frac{x^n}{e^x+e^{-x}}= \beta(n+1) \Gamma(n+1)}$$
