We have simulated Brownian bead-spring chains of up to 125 units with fluctuating hydrodynamic and excluded volume interactions using the Chebyshev polynomial approximation proposed by Fixman [Macromolecules 19, 1204 (1986)] for the square root of the diffusion tensor. We have developed a fast method to continuously determine the validity of the eigenvalue range used in the polynomial approximation, and demonstrated how this range may be quickly updated when necessary. We have also developed a weak first order semiimplicit time integration scheme which offers increased stability in the presence of steep excluded volume potentials. The full algorithm scales roughly as O(N-2.25) and offers substantial computational savings over the standard Cholesky decomposition. The above algorithm was used to obtain scaling exponents for various static and zero shear rate dynamical properties, which are found to be consistent with theoretical and/or experimental predictions.