Power flow calculations are crucial for the study of power systems, as they can be used to calculate bus voltage magnitudes and phase angles, as well as active and reactive power flows on lines. In this paper, a new approach, the Recycling Newton–Krylov (ReNK) algorithm, is proposed to solve the linear systems of equations in Newton–Raphson iterations. The proposed method uses the Generalized Conjugate Residuals with inner orthogonalization and deflated restarting (GCRO-DR) method within the Newton–Raphson algorithm and reuses the Krylov subspace information generated in previous Newton runs. We evaluate the performance of the proposed method over the traditional direct solver (LU) and iterative solvers (Generalized Minimal Residual Method (GMRES), the Biconjugate Gradient Stabilized Method (Bi-CGSTAB) and Quasi-Minimal Residual Method (QMR)) as the inner linear solver of the Newton–Raphson method. We use different test systems with a number of busses ranging from 300 to 70000 and compare the number of iterations of the inner linear solver (for iterative solvers) and the CPU times (for both direct and iterative solvers). We also test the performance of the ReNK algorithm for contingency analysis and for different load conditions to simulate optimization problems and observe possible performance gains. •A novel recycling Krylov subspace approach for the power flow algorithm•Investigation of the angles between the Krylov subspaces generated from the Jacobians•ILU preconditioner is used once with AMD reordering for iterative solutions•Computational improvement due to similar subspaces for contingencies/varying loads•Extensive set of experiments for several different iterative and direct solvers.