Exploiting symmetry in the axisymmetric Jeans equations I'm solving the axisymmetric Jeans equations to determine the initial conditions for an exponential disk inside an NFW halo. The density profiles of the two components are 
\begin{equation}
\rho_d(R,z)=\rho_{0,d}\exp\left(-\frac{R}{R_0}\right)\exp\left(-\frac{|z|}{z_0}\right)
\end{equation}

\begin{equation}
\rho_h(r)=\frac{\rho_{0,h}}{r/r_s(1+r/r_s)^2}
\end{equation}

Due to the axial symmetry, for the meridional velocity dispersion we have
\begin{equation}
\sigma^2(R,z)=\frac{1}{\rho}\int_z^\infty \rho \frac{\partial \Phi}{\partial z'}dz'
\end{equation}

As the galaxy is symmetric about the disk plane ($z=0$), we would expect the distribution function at a distance $h$ above the disk plane to be identical to that at $z=-h$; in other words, if we flip the galaxy, the only difference will be a change in direction of particle motion, but not speed or dispersion.

Therefore, can't we replace the $z$ in the limit with $|z|$?
Background: If $f(\mathbf{x},\mathbf{v},t)$ is the distribution function of the stars in phase space and $n = \int f\,\mathrm{d}^3\mathbf{v}$ is the star density, then in the cylindrical coordinates $(R,\varphi,z)$, the $z$-component of the second of the Jeans equations applied to an axisymmetric system is
$$\renewcommand{\expv}[1]{\langle{#1}\rangle}\partial_t(n\expv{v_z}) + \partial_R(n\expv{v_Rv_z}) + \partial_z(n\expv{v_z^2}) + \frac{n\expv{v_Rv_z}}{R} + n\partial_z\Phi = 0\text{.}$$
In a steady state, the first term vanishes, so if the positive and negative terms of the $z$-component of velocity terms balance, then so do the second and fourth, in which case:
$$\partial_z(n\expv{v_z^2}) = -n\partial_z\Phi\text{,}$$
which we can then integrate.

If the the star and mass distributions are axisymmetric and symmetric about the equatorial ($z = 0$) plane, then the potential has even symmetry in $z$, $\Phi(R,z) = \Phi(R,-z)$, and so does $n$. Thus their derivatives with respect to $z$ are odd functions, and $\Phi_{,z}(R,z) = -\Phi_{,z}(R,-z)$ implies
$$\frac{1}{n}\int_{-z}^z n\Phi_{,z'}\,\mathrm{d}z' = 0\text{,}$$
since the integrand has odd symmetry. Hence,
$$\frac{1}{n}\int_z^\infty n\Phi_{,z'}\,\mathrm{d}z' = \frac{1}{n}\int_{|z|}^\infty n\Phi_{,z'}\,\mathrm{d}z'\text{.}$$