The increasing penetration of Distributed Generation (DG) has drastically changed the topology of distribution systems. In order to fully consider this new scenario in steady-state studies, there is a need for an accurate load flow method that (a) models three phase unbalanced feeders and loads appropriately, and (b) uses a Jacobian so that the state of the system can be ascertained using analytical matrix methods. Many modern load flow techniques have been proposed for solving distribution systems. However, the majority of these methods fail in providing a clear insight on the network states due to the absence of Jacobian matrix computation. Conversely, traditional Newton-Raphson (NR) methods provide the calculation of the Jacobian matrix, but these methods have lost popularity due to their poor convergence in distribution system studies. Using recent developments of the single-phase NR solutions, an improved three-phase NR-based load flow solution is proposed to study unbalanced distribution systems and is presented in this paper. The proposed method is used to solve the IEEE test systems consisting of 4, 13, 34, 37, and 123 buses. Convergence performance is discussed in detail.