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.

\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}$

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\\ =& \>\text{Im}\int_0^1 \frac{(-1)^n i \ln^n t}{1-i t}dt=n!\>\text{Im}\>\text{Li}_{n+1}(i) \end{align}