This study presents a mathematical and computational analysis of a fractional-order tumor-immune system that includes cancer stem cell dynamics and chemotherapeutic intervention. A four-compartment model is developed utilizing the Caputo fractional derivative to characterize the interactions between cancer stem cells (S), effector immune cells (E), malignant tumor cells (T), and chemotherapeutic agent concentration (M). The essential qualitative characteristics of the model, such as existence, uniqueness, positivity, boundedness, and Ulam–Hyers stability of solutions, are meticulously demonstrated to guarantee both mathematical and biological coherence. The tumor-free equilibrium is examined, and the basic reproduction number \((R_0)\) is established as a critical parameter influencing tumor persistence and elimination. Stability analysis proves that the tumor-free equilibrium is globally asymptotically stable when \(R_0 < 1\). Numerical simulations, performed using the second-order Fractional Adams-Bashforth-Moulton (FABM) method, demonstrate that the fractional order \(\gamma\) significantly modulates the system’s kinetic behavior. Moreover, the findings indicate that reducing the fractional order results in a pronounced damping effect and increased physiological latency, reflecting the memory inherent in cellular environments. Three-dimensional sensitivity surfaces illustrate the collaborative effect of the tumor growth rate r and the chemotherapeutic killing rate \(k_T\) as the primary drivers of therapeutic success. Overall, the findings demonstrate that the presented model provides theoretical insights into tumor-immune dynamics and may assist in understanding the influence of memory effects on treatment outcomes better than classical integer-order models, offering potential insights for the optimization of therapeutic outcomes and system stability.
