Solving an integral involving exponential functions $\int_{-\infty}^{+\infty} \frac{1}{\left(e^x+ e^{-x}\right)^n} e^{-\rho x^2 + a x} dx$ I need to calculate
$$\int_{-\infty}^{+\infty} \frac{1}{\left(e^x+ e^{-x}\right)^n} e^{-\rho x^2 + a x} dx$$ where $n \in \mathbb{N}$, $\rho > 0$ and $a \in \mathbb{R}$, but I don't know how to follow. I've tried to include the expression in symbolic software trying to get a result with respect other functions, but nothing. I start thinking about approximating the integral using numerical integration, but before, I would like to be sure that the integral can not be expressed with respect other functions. Does anybody knows if I should go directly for numerical integration?
I am really lost, thank you in advance.
Updates:

Some particular cases can be computed using Wolfram Alpha. It seems that
$$\int_{-\infty}^{+\infty} \frac{1}{\left(e^x+ e^{-x}\right)} e^{-\rho x^2 +  x} dx = \frac{\sqrt{\pi/\rho}}{2}$$
and
$$\int_{-\infty}^{+\infty} \frac{1}{\left(e^x+ e^{-x}\right)} e^{-\rho x^2 + 3 x} dx = \int_{-\infty}^{+\infty} \frac{1}{\left(e^x+ e^{-x}\right)} e^{-\rho x^2 - 3 x} dx = (2 e^{1/\rho}-1)\frac{\sqrt{\pi/\rho}}{2}.$$

More generally (thank you JanG) we have
$$\int_{-\infty}^{+\infty} \frac{1}{\left(e^x+ e^{-x}\right)} e^{-\rho x^2 + (2m+1) x} dx =  \frac{\sqrt{\pi/\rho}}{2} (-1)^m \left( 1 + 2 \sum_{\ell=1}^m (-1)^\ell e^{\ell^2/\rho}\right)$$
for $m$ non-negative integer.
 A: Integral in general is terrible, so in practice I would suggest numerical methods. However, it can be represented as a series in a complicated way.
First, let's deal with the power in the denominator. We have:
$$\int_{-\infty}^{+\infty} \frac{e^{-\rho x^2 + a x}}{\left(e^x+ e^{-x}\right)^n}  dx=\int_{-\infty}^{+\infty} \frac{e^{-\rho x^2 + (a+n) x}}{\left(1+ e^{2x}\right)^n}  dx$$
Let's introduce a new parameter $b$:
$$I_n(\rho,a,b)=\int_{-\infty}^{+\infty} \frac{e^{-\rho x^2 + (a+n) x}}{\left(b+ e^{2x}\right)^n}  dx$$
Introducing a new integral:

$$J(\rho,a,b,n)=\int_{-\infty}^{+\infty} \frac{e^{-\rho x^2 + (a+n) x}}{b+ e^{2x}}  dx$$

We have:

$$I_n=\frac{(-1)^{n-1}}{(n-1)!} \frac{\partial^{n-1} }{\partial~ b^{n-1}} J \tag{1}$$

Thus, we have simplified the problem to finding $J$ in analytic form. It's still too complicated to obtain a closed form (as far as I know), but the series solution is possilbe.
Let $y=2x$, $\rho=4R$ and $a+n=2A$:

$$J=\frac{1}{2} \int_{-\infty}^{+\infty} \frac{e^{-R y^2 + A y}}{b+ e^y}  dy=\frac{1}{2} \int_{-\infty}^{\ln b} \frac{e^{-R y^2 + A y}}{b+ e^y}  dy+\frac{1}{2} \int_{\ln b}^{+\infty} \frac{e^{-R y^2 + A y}}{b+ e^y}  dy$$

Separation of the limits allows us to represent each denominator as a series:
$$\frac{1}{b+ e^y}=\sum_{k=0}^\infty (-1)^k b^{-k-1} e^{ky}, \qquad y< \ln b$$
$$\frac{1}{b+ e^y}=\sum_{l=0}^\infty (-1)^l b^l e^{-(l+1)y}, \qquad y> \ln b$$
Skipping the standard integration we obtain:

$$J=\frac{1}{4} \sqrt{\frac{\pi }{R}} \left[\sum _{k=0}^{\infty } (-1)^k b^{-k-1} e^{\frac{(A+k)^2}{4 R}} \left(\text{erf}\left[\sqrt{R} \ln b-\frac{A+k}{2 \sqrt{R}}\right]+1\right)+ \\ + \sum _{l=0}^{\infty } (-1)^l b^l e^{\frac{(A-l-1)^2}{4 R}} \left(\text{erf}\left[\frac{A-l-1}{2 \sqrt{R}}-\sqrt{R} \ln b \right]+1\right)\right] \tag{2}$$

Here $\text{erf}$ is the error function.
Using $(1)$ and $(2)$ we obtain a (terrible) analytic expression for the general integral $I_n$.
