The azimuthal-invariant, 3-d cylindrical, incompressible Navier-Stokes equations are solved to steady state for a finite-length, physically realistic model. The numerical method relies on an alternating-direction implicit (ADI) scheme that is formally second-order accurate in space and first-order accurate in time. The equations are linearized and uncoupled by evaluating variable coefficients at the previous time iteration. Wall grid clustering is provided by a Roberts transformation in radial and axial directions. A vorticity-velocity formulation is found to be preferable to a vorticity-streamfunction approach. Subject to no-slip, Dirichlet boundary conditions, except for the inner cylinder rotation velocity (impulsive start-up) and zero-flow initial conditions, nonturbulent solutions are obtained for sub- and supercritical Reynolds numbers of 100 to 400 for a finite geometry where R[subscript] outer/R[subscript] inner = 1.5, H/R[subscript] inner = 0.73 and H/[delta]R = 1.5. An axially-stretched model solution is shown to asymptotically approach the 1-d analytic Couette solution at the cylinder midheight. Flowfield change from laminar to Taylor-vortex flow is discussed as a function of Reynolds number. Three-dimensional velocities, vorticity and streamfunction are presented via 2-d graphs and 3-d surface and contour plots. A Prandtl-Van Driest turbulence model based on an effective isotropic eddy viscosity hypothesis was applied resulting in accurate 1-d turbulent flow solutions assuming long cylinders. A small aspect ratio correction factor was empirically determined. Comparisons to experiment are very good. Extending the nonturbulent analysis, 3-d turbulent flow equations are developed for Prandtl-Van Driest and energy-dissipation turbulence models. The energy-dissipation model includes corrections for streamline curvature, system rotation and low-Re effects. Solutions of the 3-d equations involve current work in progress.