This paper develops an innovative derivative-free optimization method for solving large-scale nonlinear monotone systems on Hadamard manifolds. The proposed algorithm is built as a convex combination of two classical conjugate gradient schemes, Fletcher–Reeves and Polak–Ribiere–Polyak, thereby inheriting the good global convergence properties of the former and the practical efficiency of the latter. Unlike conventional Riemannian optimization methods, which typically require gradients or Jacobians that may be unavailable or prohibitively expensive to compute on curved spaces, the new approach is genuinely derivative-free. It integrates hybrid conjugate gradient strategies with hypersurface projection techniques defined via retractions and vector transport, ensuring that all iterates remain on the manifold while preserving suitable conjugacy and descent properties. One feature of the algorithm is its function-based line search, which employs Armijo- and Wolfe-type conditions to guarantee sufficient descent without using derivatives of the underlying operator, thereby mitigating stagnation issues commonly observed in derivative-free schemes. Under standard assumptions of monotonicity and Lipschitz continuity, a rigorous convergence analysis establishes global convergence of the iterates to a solution of the monotone system. Extensive numerical experiments on test problems with dimensions up to 50,000 illustrate the method’s efficiency, robustness, and competitiveness relative to existing state-of-the-art algorithms for large-scale nonlinear systems. Finally, when applied to an image restoration problem, the method successfully reconstructs severely degraded images, achieving high-quality reconstructions as quantified by standard performance metrics such as Peak Signal-to-Noise Ratio (PSNR) and Signal-to-Noise Ratio (SNR), demonstrating its practical relevance in imaging applications.