This paper proposes a convex optimization based method that either locates all real roots of a set of power flow equations or declares no real solution exists in the given area. In the proposed method, solving the power flow equations is reformulated as a global optimization problem (GPF for short) that minimizes the sum of slack variables. All the global minima of GPF with a zero objective value have a one-to-one correspondence to the real roots of power flow equations. By solving a relaxed version of GPF over a hypercube, if the optimal value is strictly positive, there is no solution in this area and the hypercube is discarded. Otherwise the hypercube is further divided into smaller ones. This procedure repeats recursively until all the real roots are located in small enough hypercubes through the successive refinement of the feasible region embedded in a bisection paradigm. This method is desired in a number of power system security assessment applications, for instance, the transient stability analysis as well as voltage stability analysis, where the closest unstable equilibrium and all Type-I unstable equilibrium is required, respectively. The effectiveness of the proposed method is verified by analyzing several test systems.