Does anyone knows how to average function in three dimensions $$
\begin{split}
\frac{dr}{dt}       &= -\epsilon \sin^2(\theta)z \\
\frac{d \theta}{dt} &= -1-\epsilon \cos\theta \sin\theta z\\
\frac{dz}{dt}       &= \epsilon(r^2-T)
\end{split}
$$
Does anyone know how to do function averaging by using Fourier analysis in cylindrical coordinates or have some helpful reference?
 A: I would refer to F. Verhulst, Nonlinear Differential Equations and Dynamical Systems (Springer, 1996), chapter 11, section 11.6. Quoting equation (11.37) therein, you have
\begin{align}
\dot{x} &= \epsilon\, X(\phi,x) + \mathcal{O}(\epsilon^2),\quad x\in D\subset \mathbb{R}^n, \tag{1a}\\
\dot{\phi} &= \Omega(x) + \mathcal{O}(\epsilon),\quad \phi \in S^1,\tag{1b}
\end{align}
with $X(\phi,x)$ periodic in $\phi$; in your case, $\phi = \theta$, $\Omega(x) = -1$, $x = (r,z)$ and $X(\phi,x) = (-\sin(\theta)^2 z, r^2 - T)$. Then, Theorem 11.4 (Verhulst, p.151) says that

Consider system $(1)$ with initial values $x(0) = x_0$, $\phi(0) = \phi_0$ and suppose that
a. the righthand sides are $C^1$ in $D \times S^1$; 
b.) the solution of
  $$
\dot{y} = \epsilon\, X^0(y),\quad y(0) = x_0
$$
  with
  $$
X^0(y) = \int_{S^1} X(\phi,x)\,\text{d}\phi
$$
  is contained in an interior subset of $D$ in which $\Omega(x)$ is bounded away from zero by a constant independent of $\epsilon$; 
then $x(t) - y(t) = \mathcal{O}(\epsilon)$ on the time-scale $1/\epsilon$.

In your case, this yields the system
\begin{align}
\dot{r}_\text{av} &= -\frac{1}{2} z_\text{av},\\
\dot{z}_\text{av} &= r_\text{av}^2 - T,
\end{align}
for the averaged values of $r$ and $z$, which is Hamiltonian and hence integrable.
