An efficient implicit Navier-Stokes method for computing steady, three- dimensional flowfields characteristic of high-speed propulsion systems is presented. A nonlinear iteration strategy based on planar Gauss-Seidel sweeps is used to drive the solution toward a steady state, with approximate factorization errors within a crossflow plane reduced by the application of a quasi-Newton technique. A hybrid discretization approach is employed, with flux-vector splitting used in the streamwise direction, and central differences with artificial dissipation used for the transverse fluxes. Convergence histories and comparisons with experimental data are presented for several three-dimensional shock/boundary-layer interactions. Turbulent closure is provided by a modification of the Baldwin-Barth one-equation model.