The nonlinear saturation of the $r$-mode instability and its effects on the spin evolution of low mass x-ray binaries (LMXBs) are modeled using the triplet of modes at the lowest parametric instability threshold. We solve numerically the coupled equations for the three modes in conjunction with the spin and temperature evolution equations. We observe that very quickly the mode amplitudes settle into quasistationary states that change slowly as the temperature and spin of the star evolve. Once these states are reached, the mode amplitudes can be found algebraically and the system of equations is reduced from eight to two equations: spin and temperature evolution. The evolution of the neutron star angular velocity and temperature follow easily calculated trajectories along these sequences of quasistationary states. The outcome depends on whether or not the star will reach thermal equilibrium, where the viscous heating by the three modes is equal to the neutrino cooling ($H=C$ curve). If, when the $r$-mode becomes unstable, the star spins at a frequency below the maximum of the $H=C$ curve, then it will reach a state of thermal equilibrium. It can then either (1) undergo a cyclic evolution with a small cycle size with a frequency change of at most 10%, (2) evolve toward a full equilibrium state in which the accretion torque balances the gravitational radiation emission, or (3) enter a thermogravitational runaway on a very long time scale of $\ensuremath{\approx}{10}^{6}\text{ }\text{ }\mathrm{years}$. If the star does not reach a state of thermal equilibrium, then a faster thermal runaway (time scale of $\ensuremath{\approx}100\text{ }\text{ }\mathrm{years}$) occurs and the $r$-mode amplitude increases above the second parametric instability threshold. Following this evolution requires more inertial modes to be included. The sources of damping considered are shear viscosity, hyperon bulk viscosity, and viscosity within the core-crust boundary layer. We vary proprieties of the star such as the hyperon superfluid transition temperature ${T}_{c}$, the fraction of the star that is above the threshold for direct URCA reactions, and slippage factor, and map the different scenarios we obtain to ranges of these parameters. We focus on ${T}_{c}\ensuremath{\gtrsim}5\ifmmode\times\else\texttimes\fi{}{10}^{9}\text{ }\text{ }\mathrm{K}$ where nonlinear effects are important. Wagoner [R. Wagoner, Astrophys. J. 578, L63 (2002).] has shown that a very low $r$-mode amplitude arises at smaller ${T}_{c}$. For all our bounded evolutions the $r$-mode amplitude remains small $\ensuremath{\sim}{10}^{\ensuremath{-}5}$. The spin frequency of accreting neutron stars is limited by boundary layer viscosity to ${\ensuremath{\nu}}_{\mathrm{max}}\ensuremath{\approx}800\text{ }\text{ }\mathrm{Hz}[{S}_{\mathrm{ns}}/({M}_{1.4}{R}_{6}){]}^{4/11}{T}_{8}^{\ensuremath{-}2/11}$. Fast rotators are allowed for $[{S}_{\mathrm{ns}}/({M}_{1.4}{R}_{6}){]}^{4/11}{T}_{8}^{\ensuremath{-}2/11}\ensuremath{\sim}1$ and we find that in this case the $r$-mode instability would be active for about 1 in 1000 LMXBs and that only the gravitational waves from LMXBs in the local group of galaxies could be detected by advanced LIGO interferometers.