The azimuthal-invariant, three-dimensional 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 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(outer)/R(inner) = 1.5, H/R(inner) = 0.73, and H/Delta-R = 1.5. An axially-stretched model solution is shown to asymptotically approach the one-dimensional 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 two-dimensional graphs and three-dimensional surface and contour plots.