  We discuss a numerical method to compute the homogeneous solutions of the
Teukolsky equation which is the basic equation of the black hole perturbation
method. We use the formalism developed by Mano, Suzuki and Takasugi, in which
the homogeneous solutions of the radial Teukolsky equation are expressed in
terms of two kinds of series of special functions, and the formulas for the
asymptotic amplitudes are derived explicitly.Although the application of this
method was previously limited to the analytical evaluation of the homogeneous
solutions, we find that it is also useful for numerical computation. We also
find that so-called "renormalized angular momentum parameter", $\nu$, can be
found only in the limited region of $\omega$ for each $l,m$ if we assume $\nu$
is real (here, $\omega$ is the angular frequency, and $l$ and $m$ are degree
and order of the spin-weighted spheroidal harmonics respectively). We also
compute the flux of the gravitational waves induced by a compact star in a
circular orbit on the equatorial plane around a rotating black hole. We find
that the relative error of the energy flux is about $10^{-14}$ which is much
smaller than the one obtained by usual numerical integration methods.
