Easier approach to $\int_0^{\infty} \frac{\mathrm{e}^{-x} \cosh(2x/5)}{1 + \mathrm{e}^{-2x}} \, \mathrm{d}x$?

Starting from @Luis Sierra's answer

$$\begin{equation} I=\frac{1}{2}\int\limits_{0}^{1} \frac{t^{\frac{2}{5}}}{1+t^{2}} \,dt +\frac{1}{2}\int\limits_{0}^{1} \frac{t^{-\frac{2}{5}}}{1+t^{2}}\,dt \end{equation}$$

Using the quite standard

$$J_a=\int_0^1 \frac {t^a}{1+t^2}\, dt=\frac{1}{4} \left(\psi \left(\frac{a+3}{4}\right)-\psi \left(\frac{a+1}{4}\right)\right)\qquad \text{if} \qquad \Re(a)>-1$$ So, rearranging, $$8I=\Big[\psi\left(\frac{17}{20}\right)-\psi\left(\frac{3}{20}\right)\Big]+\Big[\psi \left(\frac{13}{20}\right)-\psi\left(\frac{7}{20}\right)\Big]=\pi \cot \left(\frac{3 \pi }{20}\right)+\pi \tan \left(\frac{3 \pi }{20}\right)$$ that is to say $$8I=\pi\csc\left(\frac{3 \pi }{20}\right)\,\sec\left(\frac{3 \pi }{20}\right)=2 \left(\sqrt{5}-1\right)\, \pi \implies I=\frac{\sqrt{5}-1}{4} \pi$$

Right before "Our contour is", you could substitute $u = v^5; du = 5v^4 dv$, and get yourself an integrand that's $$ \frac{5v^2}{1 + v^{10}}, $$ after which all you need to do is factor a tenth-degree polynomial and do endless partial-fractions work. I mean...it is elementary, just really unpleasant.

Let

\begin{equation} I=\int\limits_{0}^{+\infty} \frac{e^{-x}\cosh\left(\frac{2x}{5}\right)}{1+e^{-2x}} \,dx \end{equation}

Now let $t=e^{-x}$, which implies that $-dt=e^{-x}\,dx$. Plugging everything in yields:

\begin{equation} I=\int\limits_{0}^{1} \frac{\cosh\left(\frac{2}{5}\ln(t)\right)}{1+t^{2}} \,dt \end{equation}

Using the exponential definition of $\cosh(x)$ and the fact that $\ln(x^{a})=a\ln(x)$, we derive that:

\begin{equation} \cosh\left(\frac{2}{5}\ln(t)\right)=\frac{t^{\frac{2}{5}}+t^{-\frac{2}{5}}}{2} \end{equation}

Thus:

\begin{equation} I=\frac{1}{2}\int\limits_{0}^{1} \frac{t^{\frac{2}{5}}}{1+t^{2}} \,dt +\frac{1}{2}\int\limits_{0}^{1} \frac{t^{-\frac{2}{5}}}{1+t^{2}}\,dt \end{equation}

Now, consider the following integral:

\begin{equation} I(a,b)=\int\limits_{0}^{1}\frac{t^{a}}{1+t^{2}}e^{-bt}\,dt \end{equation}

We can compute the two integrals above with this generalized integral, note that: \begin{equation} I=\frac{1}{2}I\left(\frac{2}{5},0\right)+\frac{1}{2}I\left(-\frac{2}{5},0\right) \end{equation}

We want to obtain a differential equation with respect to $I(a,b)$, so that when we solve it, we can compute $I$. Let's take the first derivative with respect to $b$:

\begin{equation} I'(a,b)=\int\limits_{0}^{1}\frac{\partial}{\partial b}\left[\frac{t^{a}}{1+t^{2}}e^{-bt}\right]\,dt=\int\limits_{0}^{1}\frac{(-t)t^{a}e^{-bt}}{1+t^{2}}\,dt \end{equation}

Let's differentiate once again:

\begin{equation} I''(a,b)=\int\limits_{0}^{1}\frac{\partial}{\partial b}\left[\frac{(-t)t^{a}e^{-bt}}{1+t^{2}}\right]\,dt=\int\limits_{0}^{1}\frac{t^{2}t^{a}e^{-bt}}{1+t^{2}}\,dt \end{equation}

If we add and substract $1$ in the $t^{2}$, we can simplify things:

\begin{equation} I''(a,b)=\int\limits_{0}^{1}\frac{(t^{2}+1-1)t^{a}e^{-bt}}{1+t^{2}}\,dt=\int\limits_{0}^{1}\frac{(t^{2}+1)t^{a}e^{-bt}}{1+t^{2}}\,dt-\underbrace{\int\limits_{0}^{1}\frac{t^{a}e^{-bt}}{1+t^{2}}\,dt}_{I(a,b)} \end{equation}

Note that the second integral is just our original $I(a,b)$. The first integral can be expressed in terms of the lower incomplete gamma function:

\begin{equation} \int\limits_{0}^{1}\frac{(t^{2}+1)t^{a}e^{-bt}}{1+t^{2}}\,dt=\int\limits_{0}^{1}t^{a}e^{-bt}\,dt \end{equation}

If we let $z=bt$, calculate and plug everything in, one gets the following integral:

\begin{equation} \frac{1}{b}\int\limits_{0}^{b}\left(\frac{z}{b}\right)^{a}e^{-z}\,dz=\frac{1}{b^{1+a}}\underbrace{\int\limits_{0}^{b}z^{a}e^{-z}\,dz}_{\gamma\left(1+a,b\right)}=\frac{1}{b^{1+a}}\gamma\left(1+a,b\right) \end{equation}

Finally, we obtain the following differential equation:

\begin{equation} I''(a,b)+I(a,b)-\frac{1}{b^{1+a}}\gamma\left(1+a,b\right)=0 \end{equation}

Given that we have differentiated $I$ exclusively with respect to $b$, then we only need to consider the dependence with respect to $b$. Also, it is known that $\gamma\left(a+1,b\right)=a\gamma(a,b)-b^{a}e^{-b}$, then:

\begin{equation} I''(b)+I(b)-\frac{1}{b^{1+a}}\left[a\gamma(a,b)-b^{a}e^{-b}\right]=0 \end{equation}

Solving the differential equation will allow us to compute $I(a,b)$, and once this is calculated, we can just plug the necessary values to determine $I$.