Currently, several methods applied in distribution systems need power flow solutions, such as evolutionary algorithms, exhaustive search, neural network training, continuation power flow, quasi-static time-series, among others. A properly modeled distribution system requires detailing and expects a three-phase network with neutrals and groundings, which help increase simulation times. For a thorough modeling, but with reduced simulation time, distribution system equivalency can be determined using a multi-area aggregation method in three-phase four-wire multi-grounded feeders, where even the modeling of neutral cables and groundings are allowed. The proposed method enables the user to represent complex network components as a reduced size subnet, achieving a considerable computational time reduction. The proposed method can be inserted into OpenDSS , or multiphase methods, using Newton - Raphson and Backward-Forward Sweep algorithms. Results show that the method is capable of mathematically finding the correct solution for all loading levels. The exactness of the proposition is evaluated using the simple NEV test feeder, IEEE 123, the IEEE 8500 Node Test Feeder, and a meshed network having unbalanced feeders with mutual coupling between different voltage levels.