In a differential equation I am solving, the final result had an integral of $$\int\frac{dx}{e^x-x}.$$ This is non-elementary, and I was wondering if I could get rid of the integral by using some special functions.
Is there any hope in doing this? It seems that Wolfram Alpha isn't able to do it.
For an approximation of$$I=\int\frac{dx}{e^x-x}$$ a series solution would be dangerous. However, you can approximate the integrand using the $[n,n]$ corresponding Padé approximant $P_n(x)$. This would write $$\frac{1}{e^x-x}=\frac {1+\sum_{k=1}^n a_n\,x^n }{1+\sum_{k=1}^n b_n\,x^n }$$ Computing the roots and factoring $$\frac{1}{e^x-x}=\frac{a_n}{b_n}\,\frac {\prod_{k=1}^n (x-r_k)}{\prod_{k=1}^n (x-s_k)}$$ Using partial fraction decomposition $$\frac{1}{e^x-x}=\sum_{k=1}^n \frac {c_k}{x-s_k}$$ $$I=\int\frac{dx}{e^x-x}=\sum_{k=1}^n c_k \log(x-s_k)$$ will, as usual, lead to some arctangents and logarithms.
To show the quality of the approximation over a limited range, consider the norms $$\Phi_n=\int_0^{2e} \left(\frac{1}{e^x-x}-P_n\right)^2\,dx$$ $$\left( \begin{array}{cc} n & \Phi_n \\ 2 & 1.06864\times 10^{-4} \\ 3 & 8.82719\times 10^{-5} \\ 4 & 8.62984\times 10^{-8} \\ 5 & 2.08495\times 10^{-9} \\ 6 & 1.11514\times 10^{-12} \\ \end{array} \right)$$
Doing the work for $n=2$ $$P_2(x)=\frac {1-\frac{1}{3}x+\frac{1}{36}x^2 } {1-\frac{1}{3}x+\frac{19}{36}x^2 }=\frac {(x-6)^2}{19 x^2-12 x+36 }$$ $$P_2(x)=\frac 1{19}\left(1+\frac{a^2-12 a+36}{(a-b) (x-a)} -\frac{b^2-12 b+36}{(a-b) (x-b)}\right)$$ where $$a=\frac{6}{19} \left(1-3 i \sqrt{2}\right)\qquad \text{and} \qquad b=\frac{6}{19} \left(1+3 i \sqrt{2}\right)$$ So $$\int P_2(x)\,dx=\frac{x-6}{19}+\frac{306}{361} \sqrt{2} \tan ^{-1}\left(\frac{19 x-6}{18 \sqrt{2}}\right)-$$ $$\frac{108}{361} \log \left(19 x^2-12 x+36\right)$$ Computing $\int_0^t$, a few numbers $$\left( \begin{array}{ccc} t & \text{approximation} & \text{"exact"}\\ 0.5 & 0.47817 & 0.47817 \\ 1.0 & 0.84296 & 0.84308 \\ 1.5 & 1.06760 & 1.06834 \\ 2.0 & 1.19279 & 1.19499 \\ 2.5 & 1.26069 & 1.26511 \\ 3.0 & 1.29726 & 1.30445 \\ \end{array} \right)$$
Edit
Suppose that you need the integral around a given point $x=a$. Doing the same, we would have $$P_2(x)=\frac{\alpha_0+\alpha_1(x-a) +\alpha_2(x-a)^2} { 1+\beta_1(x-a) +\beta_2(x-a)^2}$$ $$\alpha_0=\frac{1}{e^a-a}\quad \alpha_1=-\frac{e^a+1}{2 \left(e^a+2\right) \left(e^a-a\right)}\quad \alpha_2=\frac{e^a}{12 \left(e^a+2\right) \left(e^a-a\right)}$$ $$\beta_1=\frac{e^a a+a+e^a+e^{2 a}-4}{2 \left(e^a+2\right) \left(e^a-a\right)}\quad \beta_2=\frac{e^a \left(-a+e^a+12\right)+6}{12 \left(e^a+2\right) \left(e^a-a\right)}$$ and the integral would be easy to compute using the same procedure to obtain (with $t=x-a$)
$$\frac{\alpha_2 }{\beta_2}t+\frac{ \left(2 \alpha_0 \beta_2^2-\alpha_1 \beta_1 \beta_2+\alpha_2 \beta_1^2-2 \alpha_2 \beta_2\right)}{\beta_2^2 \sqrt{4 \beta_2-\beta_1^2}}\tan ^{-1}\left(\frac{\beta_1+2 \beta_2 t}{\sqrt{4 \beta_2-\beta_1^2}}\right)+$$ $$\frac{(\alpha_1 \beta_2-\alpha_2 \beta_1)}{2 \beta_2^2} \log \left(\beta_1 t+\beta_2 t^2+1\right)$$
For example, using $a=\pi$ and integrating between $\frac \pi 2$ and $\frac {3\pi}2$, the above will give $0.257755$ while the solution is $0.258772$; this corresponds to a relative error of $0.39$%.
Be sure that we could do much better.