For notational convenience, let me write
[tex]
\int \frac{1}{\sqrt{2E-2(\frac{1}{r}-\frac{1}{2})+e^{-r^2-z^2}}} dz
=
\int \frac{1}{\sqrt{c + e^{-r^2} e^{-z^2}}} \, dz
[/tex]
If you plot the graph, you will see that it gets close to its asymptote y = 1/sqrt(c) rather quickly. So you could approximate it by [itex]1/\sqrt{c}[/itex], and expand around z = 0 in a series of which you only include the first several terms (they are polynomials in z with complicated coefficients, depending just on c and r though). I.e. something like
[tex]\int_{-a}^a \frac{1}{\sqrt{c + e^{-r^2} e^{-z^2}}} \, dz
=
\int_{-a}^{-\delta} \frac{1}{\sqrt{c + e^{-r^2} e^{-z^2}}} \, dz
+
\int_{-\delta}^{\delta} \frac{1}{\sqrt{c + e^{-r^2} e^{-z^2}}} \, dz
+
\int_{\delta}^a \frac{1}{\sqrt{c + e^{-r^2} e^{-z^2}}} \, dz
[/tex]
where you choose [itex]\delta \approx 2[/itex] (for sufficiently large a >> delta) conveniently so that the integrand is as good as constant outside the interval [itex][-\delta, \delta][/itex].
Then
[tex]\int_{-a}^a \frac{1}{\sqrt{c + e^{-r^2} e^{-z^2}}} \, dz
\approx
\frac{1}{\sqrt{c}} \cdot (a - \delta)
+
\int_{-\delta}^{\delta} \left( c_0 + c_2 z^2 + c_4 z^4 + \cdots \right) \, dz
+
\frac{1}{\sqrt{c}} \cdot (a - \delta)
[/tex]
where the c_{2i} are the Taylor series coefficients, e.g.
[tex]c_0 = \frac{1}{\sqrt{c + e^{-r^2}}}[/tex];
[tex]c_2 = \frac{e^{-r^2}}{2 \left(c+e^{-r^2}\right)^{3/2}}[/tex];
[tex]c_4 = \frac{1-2 c e^{r^2}}{8 \sqrt{c+e^{-r^2}} \left(c e^{r^2}+1\right)^2}[/tex]
etc. (you can get it as accurate as you want by including more terms).
Note that for [itex]a \to \infty[/itex] the integral will diverge, as the leading contribution is something like
[tex]\frac{2a}{\sqrt{c}}[/tex]