This study presents a novel fifth-order iterative method for solving nonlinear systems derived from a modified combination of Jarratt and Newton schemes, incorporating a frozen derivative of the Jacobian. The method is applied to approximate solutions of the nonlinear convection–diffusion equation. A MATLAB script function was developed to implement the approach in two stages: first, discretizing the equation using the Crank–Nicolson Method, and second, solving the resulting nonlinear systems using Newton’s iterative method enhanced by a three-step Jarratt variant. A comprehensive analysis of the results highlights the method’s convergence and accuracy, comparing the numerical solution with the exact solution derived from linear parabolic partial differential transformations. This innovative fifth-order method provides an efficient numerical solution to the nonlinear convection–diffusion equation, addressing the problem through a systematic methodology that combines discretization and nonlinear equation solving. The study underscores the importance of advanced numerical techniques in tackling complex problems in physics and mathematics.

Nonlinear systems of equations arise in a wide range of applications in engineering, physics, and biological modeling; however, classical Newton-type methods may fail when the initial approximation lies outside the basin of attraction of the desired solution. This work proposes a two-stage hybrid framework that couples Particle Swarm Optimization (PSO) for global exploration with the fifth-order Newton-Jarratt (NJN) iterative method for local refinement. The fifth-order convergence of the NJN phase, established through a complete Fréchet-derivative Taylor expansion with explicitly computed error constants, guarantees rapid local convergence once PSO delivers a sufficiently close starting point. The framework is validated on four test problems of increasing dimension (n=2,5,20,40): a two-dimensional benchmark algebraic system, a five-dimensional metabolic network model for ethanol production in Saccharomyces cerevisiae, and two large-scale systems arising from the discretization of a Hammerstein nonlinear integral equation. Over 30 independent runs per method, PSO-NJN achieves a 100% convergence rate in all four problems, with mean final residuals on the order 10⁻¹⁴-10⁻¹⁶. In comparison, pure PSO fails completely (0% success) on the high-dimensional Hammerstein cases (n=20,40) and achieves only 10% success on the metabolic model. These results confirm that combining global metaheuristic search with high-order local refinement yields a robust, scalable solver suitable for complex biological and engineering nonlinear systems.

Nonlinear systems of equations arise across engineering, physics, and biological modeling; however, classical Newton-type methods may fail when the initial approximation lies outside the convergence region of the NJN local solver. This work proposes a two-stage hybrid framework that couples Particle Swarm Optimization (PSO) for global exploration with the fifth-order Newton–Jarratt (NJN) iterative method for local refinement. The fifth-order convergence of the NJN phase, established through a complete Fréchet-derivative Taylor expansion with explicitly computed error constants, guarantees rapid local convergence once PSO delivers a sufficiently close starting point. The framework is validated on four test problems with increasing numbers of dimensions ($n=2,5,20,40$): a two-dimensional benchmark algebraic system, a five-dimensional metabolic network model for ethanol production in Saccharomyces cerevisiae, and two large-scale Hammerstein integral equation systems. Over 30 independent runs per method and under the tested conditions, PSO-NJN achieves 100% convergence with mean final residuals of order ${10}^{-14}$–${10}^{-16}$, while pure PSO fails completely on the high-dimensional Hammerstein cases ($n=20,40$) and achieves only 10% success on the metabolic model. These results confirm that combining global metaheuristic search with high-order local refinement yields a robust, scalable solver for complex biological and engineering nonlinear systems, though performance on problems with dense high-dimensional Jacobians may require further adaptation.

We develop a Laurent-series framework for the analytical characterization of singularities in any three-step frozen-Jacobian iterative operator parametrized by $(\lambda ,\alpha ,\beta )$ (where $\lambda$ controls the predictor step and $(\alpha ,\beta )$ are the weight-function coefficients) and apply it to the fifth-order Newton–Jarratt–Newton (NJN) method. Three principal results are established. First, any operator with a quadratic weight function ($\beta \ne 0$) possesses a pole of order $m=d(d+1)+1$ at every simple critical point of f, where d is the asymptotic polynomial degree, independently of $\lambda$ and $\alpha$; when $\beta =0$ the order reduces to $m=d+1$, recovering the classical Jarratt-class result. Second, for each fixed $\lambda$ there exists a unique parameter pair $({\alpha }^{*}\left(\lambda \right),{\beta }^{*}\left(\lambda \right))$ that minimises the Laurent residue $|a_{-1}|$ subject to fifth-order convergence; for the NJN ($\lambda =2/3$) the explicit solution ${\alpha }^{*}=\frac{1}{4}$, ${\beta }^{*}=\frac{3}{8}$ provides the first analytical justification for these design values. Third, the Laurent Validity Radius ${\rho }_L={\inf }_{z_s\in S\setminus \left\{z_c\right\}}|z_c-z_s|$ is the exact boundary of the domain in which the single-pole representation is valid; starting points satisfying $|z_0-z_c|\ll \left|\kappa \right|$ experience pole-dominated first-iterate growth $|G\left(z_0\right){|\sim |\kappa |}^m{\delta }^{-m}\to +\infty$. All results are validated on five test functions (including polynomial, transcendental, meromorphic, and non-symmetric cases) and all three, quantities can be computed from $f\left(z_c\right)$, $f^{\prime }\left(z_c\right)$, $f^{″}\left(z_c\right)$ alone.