There are different methodologies applied to simulate the inrush phenomenon in power transformers. Generally, these differ from each other with regard to the technique used to model the saturation curve, the numerical method used to discretize the system and the technique used to solve the problem. Depending on the techniques and methods employed in these processes, the solution to the problem may or may not present problems regarding accuracy, numerical stability, and computational efficiency. In this work, a methodology is proposed so that it eliminates or at least reduces the three previously mentioned problems. This methodology is described step-by-step. The proposal consists in modeling the saturation curve by true nonlinear equations, discretizing the system using the second-order BDF method and obtaining the solution using the Newton-Raphson method. This proposal proved to be effective with regard to the three parameters mentioned. In addition, examples of limitations arising from known methodologies were presented, in order to elucidate the advantages of the proposal.