  The manuscript addresses the problem of finding all solutions of power flow
equations or other similar nonlinear system of algebraic equations. This
problem arises naturally in a number of power systems contexts, most
importantly in the context of direct methods for transient stability analysis
and voltage stability assessment. We introduce a novel form of homotopy
continuation method called the numerical polynomial homotopy continuation
(NPHC) method that is mathematically guaranteed to find all the solutions
without ever encountering a bifurcation. The method is based on embedding the
real form of power flow equation in complex space, and tracking the generally
unphysical solutions with complex values of real and imaginary parts of the
voltage. The solutions converge to physical real form in the end of the
homotopy. The so-called $\gamma$-trick mathematically rigorously ensures that
all the paths are well-behaved along the paths, so unlike other continuation
approaches, no special handling of bifurcations is necessary. The method is
\textit{embarrassingly parallelizable} and can be applied to reasonably large
sized systems. We demonstrate the technique by analysis of several standard
test cases up to the 14-bus system size. Finally, we discuss possible
strategies for scaling the method to large size systems, and propose several
applications for transient stability analysis and voltage stability assessment.
