A High-Order Parallel Framework for Simultaneous Root-Finding in Nonlinear Systems with Multiple Solutions

Nonlinear systems with multiple roots arise frequently in biomedical and engineering models, yet their reliable numerical solution remains a challenging task. Many classical methods suffer from sensitivity to initial guesses, reduced convergence rates, and loss of accuracy in the presence of multiple or clustered solutions. In addition, the exploitation of parallelism to improve robustness and computational efficiency has received limited attention. In this work, we propose a high-accuracy parallel numerical framework of fourth-order convergence for the simultaneous approximation of all solutions of nonlinear systems with multiple roots. The proposed scheme is derivative-free and structurally decoupled, enabling efficient parallel implementation and robust convergence even when reliable initial approximations are unavailable. The effectiveness of the method is demonstrated on representative biomedical engineering models, including a glucose–insulin–glucagon regulatory network and a multi-compartment pharmacokinetic system, both exhibiting strong nonlinearity and multistability. Numerical experiments confirm stable convergence toward distinct solution clusters, machine-level accuracy, reduced residual norms, and improved computational performance when compared with existing approaches. These results indicate that the proposed framework provides a reliable and efficient alternative for solving nonlinear systems with multiple roots in complex applied settings.