  Techniques are developed here for evaluating the r-modes of rotating neutron
stars through second order in the angular velocity of the star. Second-order
corrections to the frequencies and eigenfunctions for these modes are evaluated
for neutron star models. The second-order eigenfunctions for these modes are
determined here by solving an unusual inhomogeneous hyperbolic boundary-value
problem. The numerical techniques developed to solve this unusual problem are
somewhat non-standard and may well be of interest beyond the particular
application here. The bulk-viscosity coupling to the r-modes, which appears
first at second order, is evaluated. The bulk-viscosity timescales are found
here to be longer than previous estimates for normal neutron stars, but shorter
than previous estimates for strange stars. These new timescales do not
substantially affect the current picture of the gravitational radiation driven
instability of the r-modes either for neutron stars or for strange stars.
