The goal of this research is to produce a robust, parabolized Navier-Stokes (PNS) code that will significantly reduce the computer time required to calculate flows about complex vehicles with embedded subsonic/separated regions. The major drawback of "current day" PNS codes is that they cannot be used to compute separated regions which occur near canopies, wing-body junctures, etc. As a result, Navier-Stokes (NS) codes are often used to compute the entire flowfield despite the fact that a PNS code requires at least one order of magnitude less computer time and storage;An innovative approach has been developed to permit a PNS code to compute embedded regions that cause upstream influence. In this approach, the embedded region is automatically detected and the streamwise extent is determined prior to the computation or while the computation is in progress. The PNS equations are then solved with an iterative (IPNS) algorithm in this region to duplicate the results that would he obtained with a NS code. Once the embedded region is computed, the algorithm returns to the standard space-marching PNS mode until the next embedded region is encountered. This method has been incorporated into NASA's upwind PNS (UPS) code and validated by applying it to several 2-D test cases. These test cases include flows over compression ramps, shock-boundary-layer interactions, flows over expansion corners, and flow over a general geometry with multiple embedded regions. The results computed using this approach are in excellent agreement with NS computations and experimental data;In addition, new correlation functions have been developed that accurately predict the streamwise extent of the embedded regions for all of the geometries considered. This is the first time that any correlation (theoretical or empirical) has been shown to accurately predict where the single-sweep PNS method is inaccurate for a wide range of flow conditions. These correlation functions in conjunction with the IPNS algorithm permit completely automatic computation of steady, laminar supersonic flowfields with embedded subsonic/separated regions using a space-marching code as the primary flow solver.