Numerical models are paramount in describing the complex physical world around us. They are often based on non-linear functions, which have to be solved in order to run a simulation. The first goal of this thesis is to compare the Jacobian-free Newton-Krylov method against the regular Newton-Raphson method for solving non-linear equations. When doing this, preconditioning was only used for the Newton-Raphson method. As a second goal, the optimal parameters for solving a coupled system with the Schur complement method are determined. For computational purposes, these analyses are performed using PETSc. The comparison between JFNK and Newton-Raphson were performed by solving a heat equation. The physical system which was analyzed is a one-dimensional radiating rod, with a heterogeneous thermal conductivity and Dirichlet boundary conditions. Firstly this rod was modelled as one single system. For this case, it was found that Newton-Raphson outperformed the JFNK method by a factor of 5-1700, depending on the type of preconditioning used. Secondly, the rod was modelled as a coupled system by solving two parts of the rod separately. In this case, it was found that the JFNK method outperformed the preconditioned Newton-Raphson method by a factor of 9.1± 0.3. It was concluded that the JFNK is only favorable over regular Newton-Raphson for coupled systems. To achieve the second goal, an incompressible Navier-Stokes coupled system was solved using the Schur complement method. The coupled system originated from incompressible flow in a back-step pipe, using a finite element method. For the Schur complement parameters, it was shown that using an approximation of the Schur complement offered a significant increase in efficiency. The momentum and pressure subsystems were analyzed separately. The momentum subsystem performed best with a relative tolerance of τ_r =10^-4.5, while the pressure subsystem performed best with a tolerance of 10^-2. With these parameters, the Schur complement method had the same computational time as the pressure-correction method, which was used as a benchmark. For future research, there are three main recommendations. Firstly, to use preconditioning for the JFNK method. This was not done in this thesis because of technical constraints, but it is theoretically possible. Secondly, if the JFNK method is preconditioned, it could be used for the incompressible Navier-Stokes simulation. Finally, there are combinations of parameters which were not tested for the Schur complement method. Trying out more combinations might result in finding an even more optimized method.