This paper introduces a binary search algorithm using second order polynomial fitting to efficiently determine the maximum power transfer to a non-controlled load when accounting for variations in Thévenin voltage magnitude due to non-linearity. This is used for voltage stability boundary monitoring of a power system in real time. The binary search with polynomial fitting (BSPF) is compared to a reference algorithm, which sweeps over different load levels, and a binary search and is shown to improve both runtime and accuracy of results. The assessment method can take advantage of parallelization, which together with the BSPF algorithm makes it possible to determine a margin for each of the 2.000 non-controlled loads in a 3.000 bus test system in less than 6 seconds. This enables early detection of voltage instability in highly dynamic future smart grid based power systems.