The study introduces a fractional mathematical model in the Caputo sense for hematopoietic stem cell-based therapy, utilizing generalized Bernoulli polynomials (GBPs) and operational matrices to solve a system of nonlinear equations. The significance of the study lies in the potential therapeutic applications of hematopoietic stem cells (HSCs), particularly in the context of HIV infection treatment, and the innovative use of GBPs and Lagrange multipliers in solving the fractional hematopoietic stem cells model (FHSCM). The aim of the study is to introduce an optimization algorithm for approximating the solution of the FHSCM using GBPs and Lagrange multipliers, and to provide a comprehensive exploration of the mathematical techniques employed in this context. The research methodology involves formulating operational matrices for fractional derivatives of GBPs, conducting a convergence analysis of the proposed method, and demonstrating the accuracy of the method through numerical simulations. The major conclusion is the successful introduction of GBPs in the context of the FHSCM, featuring innovative control parameters and a novel optimization technique. The study also highlights the significance of the proposed method in providing accurate solutions for the FHSCM, thus contributing to the field of mathematical modeling in biological and medical research.