With the widespread adoption of bilateral power supply modes, the operational characteristics of traction power systems have become increasingly complex, posing higher demands on voltage stability analysis. Conventional power flow calculation methods exhibit limitations in modeling heterogeneous load structures and handling varying operating conditions, leading to constrained modeling capability as well as compromised accuracy and robustness. To address these challenges, a unified modeling method for power flow calculation suitable for bilateral supply systems is proposed. This approach integrates the modeling frameworks of three-phase grid-side loads and single-phase traction-side loads, thereby overcoming the difficulties of the Newton-Raphson method in obtaining load power on the catenary and rail sides, as well as the limitations of the continuous linear method in terms of insufficient accuracy and inability to dynamically update the Jacobian matrix. On this basis, an extended continuation power flow model is developed to systematically analyze the impact of different supply modes, line structures, and short-circuit capacities on system operational characteristics. Simulation results demonstrate that the proposed method accurately generates P-V curves and identifies critical voltage stability points, offering significant advantages in model universality, computational accuracy, and engineering applicability.