Load flow (LF) analysis is one of the most important aspects in power system studies. It is the most significant and necessary way to investigate the problems in power system operation and planning. The LF problem comprises a set of nonlinear algebraic equations that must be solved mathematically through iterations. The solution convergence is the most important criteria that is largely affected by the size of power system, which continues to increase in the current modern power system field. Thus, there is no guarantee that the iterative approaches will converge to a valid solution for problems with such large dimensions. This paper develops generalized, effective and simple mathematical models for the solution of the LF problem. The problem is first solved by Gauss–Seidel (GS) method in order to generate the training data. Eureqa software is then adopted for the purpose of data training and mathematical models generation. The models relate bus voltage magnitudes and angles as output parameters with the load active and reactive power values as input parameters. To study the validity of the proposed approach, the mathematical models have been developed for two benchmark test systems; IEEE 5-bus test system and 9-bus test system developed by Western Systems Coordinating Council (WSCC). When compared with the outputs resulting from GS technique, the results have shown efficient and accurate capability of the generated models for evaluating bus voltages magnitudes and angles as well as generated reactive power. The models have also been compared with other published research. The results have shown efficient performance of the developed models.