This function solves a system of non-linear equations using the Jacobian-Free Newton-Krylov (JFNK) method. The main advantage of using JFNK over the traditional Newton method is to avoid the need for generating and inverting the Jacobian matrix. Typically the Jacobian matrix is not analytically attainable and its numerical approximation (e.g., via finite-difference methods) is not easily invertible.

Reference: (1) Knoll D.A., Keyes D.E. "Jacobian-Free Newton-Krylov methods: a survey of approaches and applications", Journal of Computational Physics, August 2003.