This dissertation centers on two major aspects dictating the computational time of applications based on the solution of systems of coupled nonlinear parabolic equations: nonlinear and linear iterations. The former aspect leads to the conception of a novel way of reusing the Krylov information generated by GMRES for solving linear systems arising within a Newton method. The approach stems from theory recently developed on a nonlinear version of the Eirola-Nevanlinna, algorithm (originally for solving non-symmetric linear systems) which is capable of converging twice as fast as Broyden's method. A secant update strategy of the Hessenberg matrix resulting from the Arnoldi process in GMRES amounts to reflecting a secant update of the current Jacobian with the rank-one term projected onto the generated Krylov subspace (Krylov-Broyden update). This allows the design of a new nonlinear Krylov-Eirola-Nevanlinna (KEN) algorithm and a higher-order version of Newton's method (HOKN) as well. The underlying development is also auspicious to replace the use of GMRES by cheaper Richardson iterations for the sake of fulfilling the inexact Newton condition. Hence, three algorithms derived from Newton's method, Broyden's method and the nonlinear Eirola-Nevanlinna algorithm are proposed as a part of a new family of hybrid Krylov-secant methods. All five algorithms are shown to be computationally more economical than their Newton and quasi-Newton counterparts. The aspect of linear iterations complements the present research with an analysis on nested or inner-outer iterations to efficiently precondition Krylov subspace iterative solvers for linear systems arising from systems of coupled nonlinear equations. These preconditioners are called two-stage preconditioners and are developed on the basis of a simple but effective decoupling strategy. Their analysis is restricted to the particular class of problems arising in multi-phase flow phenomena modeled by systems of coupled nonlinear parabolic equations. The resulting approach outperforms fairly robust and standard preconditioners that "blindly" precondition the entire coupled linear system. Theoretical discussion and computational experiments show the suitability that both linear and nonlinear aspects undertaken in this research have for large scale implementations.