EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6166 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Investigation Based on the Chebyshev-HPM for Breast Cancer as a Mathematical Model M. Adel1, M. M. Khader2,∗, Mohammed Messaoudi2 1 Department of Mathematics, Faculty of Science, Islamic University of Madinah, Medina, Saudi Arabia 2 Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia Abstract. We present the approximate solution for mathematical model for Breast Cancer (BC) over three time intervals. The suggested approach depends on the homotopy perturbation method developed with Chebyshev series (CHPM). The residual error function is calculated and used as a basic criterion in evaluating the efficiency/accuracy of the presented numerical scheme. We utilize the solution by RK4 method for comparison with the results of the method used. Through these results, we can confirm that the applied technique is an effective tool to give a simulation of such models. Illustrative instances are given to confirm the validity and usefulness of the proposed procedure. 2020 Mathematics Subject Classifications: 34A08, 65N06 Key Words and Phrases: Breast cancer, homotopy perturbation method, Chebyshev expansion, RK4 method 1. Introduction Breast cancer is characterized as a pathological disorder arising from the unregulated growth of cells in breast tissue. According to extensive statistics from the WHO regarding the worldwide burden of cancer, BC exhibits the highest spread rate relative to other types of cancer [1]. According to studies done by the WHO, BC was classified as the sec- ond most widespread cancer in 2004, posing a significant threat to women since it affects around 8-9 percent of the global female population [2]. Notwithstanding extensive research and numerous inquiries, the precise etiology of BC remains ambiguous ([2], [3]). BC is a widespread malignancy among women post-puberty, with its prevalence escalating with advancing age [2]. In light of the factors above, we underscore the necessity for a thor- ough comprehension of the epidemiology of BC and its implications for women’s health, ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6166 Email addresses: adel@sci.cu.edu.eg (M. Adel), mmkhader@imamu.edu.sa (M. M. Khader), mmessaoudi@imamu.edu.sa (M. Messaoudi) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 2 of 16 since this information is critical for formulating efficacious precautionary and therapeutic strategies globally [4]. Mathematical modeling is crucial for comprehending and investigating cancer tumors broadly, with a specific emphasis on breast cancer in this research paper. It serves to characterize and simulate tumor growth and behavior, as well as their interactions with adjacent tissues and the immune system ([5], [6]). Researchers and doctors can use these models to better understand how tumors grow, predict how treatments will work, and come up with better ways to use treatments ([7]-[9]). For a more concentrated examination of such systems, including their applications and properties, see ([10]-[13]). Due to the problems of accuracy and convergence faced by most numerical methods, a large number of researchers have used the semi-analytical method (SAM) to investigate such problems and gain deeper insights into their complexities. The advantage of the SAM as the preferred method for finding analytical approximations to complex problems contain highly nonlinear terms ([14], [15]) lies in the challenges associated with obtaining exact solutions using traditional analytical techniques. Among these methods, which have attracted the attention of many scientists and have found application in solving a wide range of complex problems, is the HPM in particular. An improved version of this method called the new HPM, was presented in 2010 in the research paper [16] and was then applied to obtain approximate solutions to the quadratic RDE. This approach has become a popular choice in many studies ([17], [18]) due to its simplicity and ease of mathematical calculation. This distinction and simplicity can be achieved by defining the first approximation as a power series and then setting all iterations to zero except for the initial iteration. That is, the power series takes a central place in the method, which in turn facilitates obtaining approximate solutions to nonlinear equations in the form of Taylor expansion. Based on what was mentioned above about the importance of using power series, whether Taylor expansion or others, and in being more effective in improving the accuracy and convergence of analytical approximation solutions, the Chebyshev series, which depends on orthogonal Chebyshev polynomials, were used to benefit from the properties of these functions in a good approximation of functions, as well as their distinction in the faster convergence rate than other functions such as Taylor series, for example, but not limited to [19]. Due to these important properties of handling nonlinear equations and providing very accurate approximations, the transformed Chebyshev series has been widely applied in many academic studies, we mention, for example, but not limited to point equations of motion in their linear/nonlinear form [20], Cauchy problems which are represented by second-order ordinary differential equations [21], nonlinear integro-differential equations of fractional order [22], space-variable approximation of the Burger equation [23], the initial and BVPs associated with the fractional heat equation [24]. The present research focuses on developing semi-analytical methods to overcome the challenges inherent in these methods, taking advantage of the pivotal role played by the Chebyshev series, as shown in the historical context mentioned above, as well as addressing the time-consuming nature of numerical methods, mitigating these challenges. The major objective of this manuscript is to apply the CHPM by combining the ef- M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 3 of 16 fective Chebyshev series (CS) with the new homotopy perturbation approach, to use a new, accurate, and efficient analytical approximation approach used for the first time to address the present problem. This new method was compared with some existing methods to confirm its efficiency and high accuracy, and perfect agreement was obtained. The paper is organized as follows: In Section 2, we give the formulation of the proposed model. Section 3 presents the procedure solution by giving the basic concepts of the Chebyshev-HPM. Section 4 gives the numerical implementation of the proposed technique. Section 5 introduces the numerical simulation of the proposed model. Finally, Section 5 gives the conclusions and remarks. 2. Mathematical model of BC The mathematical modeling of biological events, such as BC, is a key tool for under- standing how tumors change during treatment and for dealing with epidemiological issues. This is evident from the studies conducted before ([25]-[27]). We investigate this specific model of the BC ([28], [29]): ψ̇1(t) = ψ1(t)λ1 (p1 − β1ψ1(t)− α1 ψ2(t))− (1− p)θ1 ψ1(t)ψ4(t), ψ1(0) = ψ̂0 1, ψ̇2(t) = ψ2(t)λ2(p2d− β2ψ2(t)− α2ψ3(t))− κψ2(t) + (1− p)θ1ψ1(t)ψ2(t)ψ4(t), ψ2(0) = ψ̂0 2, ψ̇3(t) = σρ+ ψ3(t)λ1(p3 − β3 − α3ψ2(t))− (1− p)θ2ψ3(t)ψ4(t), ψ3(0) = ψ̂0 3, ψ̇4(t) = α4ψ4(t) + ϱ(1− p), ψ4(0) = ψ̂0 4. (1) This section will explain what the parameters (∈ R+) of the model (1) mean in more detail ([26], [30]). The stability and convergence analysis of the model in question are thoroughly detailed in [29]. 3. Basic concepts of the CHPM The novel approach mainly depends on the use of the Chebyshev series in the new HPM. Here, we mention some basic assuming of the Chebyshev series, the new homotopy perturbation algorithm. 3.1. Chebyshev series If we have a continuous function q(t) in [a, b], then it can be rewritten in terms of the CS of the first kind as follows ([31]-[33]): q(t) = ∞∑ k=0 ′ ckTk ( 2t− b− a b− a ) , (2) M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 4 of 16 where ′ sign indicates that the coefficient of T0(t) must be reduced by half, Tk(t) = cos(k cos−1(t)) where ck = 2 π ∫ 1 −1 (1− t2)−0.5 q(0.5(b+ a+ (b− a)t))Tk(t)dt. By the following recurrence relation, we can obtain Chebyshev polynomials of the first kind given on [−1, 1]: Ts(t) = 2t Ts−1(t)− Ts−2(t), s = 2, 3, ..., T0(t) = 1, T1(t) = t. (3) Furthermore, these functions can be expressed analytically using the finite sum of powers of t as follows: Tk(t) = ⌊ k 2 ⌋∑ i=0 (−1)i 2k−2i−1 k k − i ( k − i i ) tk−2i. (4) If the domain of the problem under study is [0, 1], we need to derive the so-called Chebyshev transformed polynomials Tk(t) of the first kind, using an appropriate linear transformation via the recurrence relation (3), as follows: Tk(t) = 2(2t− 1)Tk−1(t)−Tk−2(t), k = 2, 3, ..., T0(t) = 1, T1(t) = 2t− 1. (5) 3.2. The new HPM To illustrate the process of this technique, we suppose a general form of the nonlinear ODEs as follows ([17], [18]): L[Υ(z)] +R[Υ(z)] +N [Υ(z)] + s(z) = 0, z ∈ [a, b], (6) N is a nonlinear term, and s(z) is a given source term. By harnessing the principle of homotopy, we get: H(Υ, ℓ) = L[Υ(z)]−Υ∗(z) + ℓΥ∗(z) + ℓ(R[Υ(z)] +N [Υ(z)] + s(z)) = 0, (7) where Υ∗(z) is an initial solution to the proposed problem & 0 ≤ ℓ ≤ 1 is an embedding parameter. From the principle of homotopy (7), we can find the following: L[Υ(z)] = Υ∗(z)− ℓΥ∗(z)− ℓ(R[Υ(z)] +N [Υ(z)] + s(z)), (8) Υ(z) = L−1[Υ∗(z)]− ℓ L−1[Υ∗(z)]− ℓ L−1[R[Υ(z)] +N [Υ(z)] + s(z)]. (9) Where: Υ(z) = ∞∑ k=0 ℓk uk(z), Υ∗(z) = ∞∑ k=0 ck vk(z), M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 5 of 16 then, we obtain: ∞∑ k=0 ℓk uk(z) = L−1 [ ∞∑ k=0 ck vk(z) ] − ℓ L−1 [ ∞∑ k=0 ck vk(z) ] − ℓ L−1 [ R [ ∞∑ k=0 ℓk uk(z) ] +N [ ∞∑ k=0 ℓk uk(z) ] + s(z) ] . (10) Equating the coefficients of ℓk, k = 0, 1, 2, ... on both sides, we obtain: ℓ0 : u0(z) = L−1 [ ∞∑ k=0 ck vk(z) ] , ℓ1 : u1(z) = −L−1 [ ∞∑ k=0 ck vk(z) ] − L−1[R(u0) +N(u0) + s(z)], ℓ2 : u2(z) = −L−1[R(u0, u1) +N(u0, u1)], · · · , ℓj : uj(z) = −L−1[R(u0, u1, · · · , uj−1) +N(u0, u1, · · · , uj−1)], j = 3, 4, · · · . (11) The exact solution to the equation under study will be found by finding the values of the unknown coefficient ck, assuming that u1 = 0, which will take the following form: Υ(z) = u0(z) = L−1 [ ∞∑ k=0 ck vk(z) ] . (12) 3.3. The CHPM algorithm To clarify the novel method’s methodology, we shall outline the following stages for its application: Step 1: We use the same principle of homotopy defined in (8). Step 2: Take the L−1 operator to Eq.(8), and we rearrange it to get Υ(z) as defined in (9). Step 3: Assume that Υ(z) = ∞∑ k=0 ℓk uk(z), Υ∗(z) = ∞∑ k=0 ck Tk(z), then, we get ∞∑ k=0 ℓk uk(z) = L−1 [ ∞∑ k=0 ck Tk(z) ] − ℓ L−1 [ ∞∑ k=0 ck Tk(z) ] − ℓ L−1 [ R [ ∞∑ k=0 ℓk uk(z) ] +N [ ∞∑ k=0 ℓk uk(z) ] + s(z) ] . (13) M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 6 of 16 Step 4: Comparing the coefficients of ℓk, k = 0, 1, 2, .. on both sides of (13), leads to find: ℓ0 : u0(z) = L−1 [ ∞∑ k=0 ck Tk(z) ] , ℓ1 : u1(z) = −L−1 [ ∞∑ k=0 ck Tk(z) ] − L−1[R(u0) +N(u0) + s(z)], ℓ2 : u2(z) = −L−1[R(u0, u1) +N(u0, u1)], · · · , ℓj : uj(z) = −L−1[R(u0, u1, · · · , uj−1) +N(u0, u1, · · · , uj−1)], j = 3, 4, · · · . (14) Step 5: We assume that u1 = 0, and find the values of ck, k = 0, 1, 2, ..,m by equating the Chebyshev polynomials and solving the resulting system of equations. Then, the analytical solution can be formulated as follows: Υm(z) = u0(z) = L−1 [ m∑ k=0 ck Tk(z) ] . (15) 4. Numerical implementation The CHPM technique is used to solve the current problem (1). The primary phases of the new technique are given below: Step 1: Applying the homotopy property on Eq.(1), we find: ψ̇1(t) + ℓψ∗ 1(t)− ψ∗ 1(t)− ℓ [ψ1(t)λ1 (p1 − β1ψ1(t)− α1 ψ2(t))− (1− p)θ1 ψ1(t)ψ4(t)] = 0, ψ̇2(t) + ℓψ∗ 2(t)− ψ∗ 2(t)− ℓ [ ψ2(t)λ2(p2d− β2ψ2(t)− α2ψ3(t)) − κψ2(t) + (1− p)θ1ψ1(t)ψ2(t)ψ4(t) ] = 0, ψ̇3(t) + ℓψ∗ 3(t)− ψ∗ 3(t)− ℓ [σρ+ ψ3(t)λ1(p3 − β3 − α3ψ2(t))− (1− p)θ2ψ3(t)ψ4(t)] = 0, ψ̇4(t) + ℓψ∗ 4(t)− ψ∗ 4(t)− ℓ [α4ψ4(t) + ϱ(1− p)] = 0. (16) Step 2: Taking L−1 = ∫ t 0 (.)dt, for both sides of (16) yields: ψ1(t) = ψ1(0)− ℓ L−1(ψ∗ 1(t)) + L−1(ψ∗ 1(t)) + ℓ L−1 [ψ1(t)λ1 (p1 − β1ψ1(t)− α1 ψ2(t))− (1− p)θ1 ψ1(t)ψ4(t)] , ψ2(t) = ψ2(0)− ℓ L−1(ψ∗ 2(t)) + L−1(ψ∗ 2(t)) + ℓ L−1 [ ψ2(t)λ2(p2d− β2ψ2(t)− α2ψ3(t))− κψ2(t) + (1− p)θ1ψ1(t)ψ2(t)ψ4(t) ] , ψ3(t) = ψ3(0)− ℓ L−1(ψ∗ 3(t)) + L−1(ψ∗ 3(t)) + ℓ L−1 [σρ+ ψ3(t)λ1(p3 − β3 − α3ψ2(t))− (1− p)θ2ψ3(t)ψ4(t)] , ψ4(t) = ψ4(0)− ℓ L−1(ψ∗ 4(t)) + L−1(ψ∗ 4(t)) + ℓ L−1 [α4ψ4(t) + ϱ(1− p)] . (17) M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 7 of 16 Step 3: Assuming that ψj(t) = ∞∑ k=0 ℓk ψj,k(t), ψ∗ j (t) = ∞∑ k=0 cj,k Tk(t), j = 1, 2, 3, 4, then, we get: ∞∑ k=0 ℓk ψ1,k(t) = ψ1(0)− ℓ L−1 ( ∞∑ k=0 c1,k Tk(t) ) + L−1 ( ∞∑ k=0 c1,k Tk(t) ) + ℓ L−1 [ λ1 ( p1 − β1 ( ∞∑ k=0 ℓk ψ1,k(t) ) − α1 ( ∞∑ k=0 ℓk ψ2,k(t) ))( ∞∑ k=0 ℓk ψ1,k(t) ) − (1− p)θ1 ( ∞∑ k=0 ℓk ψ1,k(t) )( ∞∑ k=0 ℓk ψ4,k(t) )] , (18) ∞∑ k=0 ℓk ψ2,k(t) = ψ2(0)− ℓ L−1 ( ∞∑ k=0 c2,k Tk(t) ) + L−1 ( ∞∑ k=0 c2,k Tk(t) ) + ℓ L−1 [ λ2 ( p2d− β2 ( ∞∑ k=0 ℓk ψ2,k(t) ) − α2 ( ∞∑ k=0 ℓk ψ3,k(t) ))( ∞∑ k=0 ℓk ψ2,k(t) ) − κ ( ∞∑ k=0 ℓk ψ2,k(t) ) + (1− p)θ1 ( ∞∑ k=0 ℓk ψ1,k(t) )( ∞∑ k=0 ℓk ψ2,k(t) )( ∞∑ k=0 ℓk ψ4,k(t) )] , (19) ∞∑ k=0 ℓk ψ3,k(t) = ψ3(0)− ℓ L−1 ( ∞∑ k=0 c3,k Tk(t) ) + L−1 ( ∞∑ k=0 c3,k Tk(t) ) + ℓ L−1 [ σρ+ λ1 ( p3 − β3 − α3 ( ∞∑ k=0 ℓk ψ2,k(t) ))( ∞∑ k=0 ℓk ψ3,k(t) ) − (1− p)θ2 ( ∞∑ k=0 ℓk ψ3,k(t) )( ∞∑ k=0 ℓk ψ4,k(t) )] , (20) ∞∑ k=0 ℓk ψ4,k(t) = ψ4(0)− ℓ L−1 ( ∞∑ k=0 c4,k Tk(t) ) + L−1 ( ∞∑ k=0 c4,k Tk(t) ) + ℓ L−1 [ α4 ∞∑ k=0 ℓk ψ4,k(t) + ϱ(1− p) ] . (21) Step 4: Equaling the terms for the equations (18)-(21) which have the same powers of ℓ: ℓ0 : ψj,0(t) = ψj(0) + L−1 ( ∞∑ k=0 cj,k Tk(t) ) , M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 8 of 16 ℓ1 : ψ1,1(t) = −L−1 ( ∞∑ k=0 c1,k Tk(t) ) + L−1 [ψ1,0(t)λ1 (p1 − β1ψ1,0(t)− α1 ψ2,0(t))− (1− p)θ1 ψ1,0(t)ψ4,0(t)] , ψ2,1(t) = −L−1 ( ∞∑ k=0 c2,k Tk(t) ) + L−1 [ ψ2,0(t)λ2(p2d− β2ψ2,0(t)− α2ψ3,0(t))− κψ2,0(t) + (1− p)θ1ψ1,0(t)ψ2,0(t)ψ4,0(t) ] , ψ3,1(t) = −L−1 ( ∞∑ k=0 c3,k Tk(t) ) + L−1 [σρ+ ψ3,0(t)λ1(p3 − β3 − α3ψ2,0(t))− (1− p)θ2ψ3,0(t)ψ4,0(t)] , ψ4,1(t) = −L−1 ( ∞∑ k=0 c4,k Tk(t) ) + L−1 [α4ψ4,0(t) + ϱ(1− p)] , (22) and so on. Step 5: We find the values cj,k, j = 1, 2, 3, 4 by assuming that ψj,1(t) = 0. Therefore, the analytical approximate solution of ψj(t) becomes as follows: ψj,m(t) = ψj,0(t) = ψj(0) + L−1 ( m∑ k=0 cj,k Tk(t) ) , j = 1, 2, 3, 4. (23) It should be noted that in step 5, the relations (Derivation, Multiplication, Integra- tion) that were mentioned earlier are applied, also, in step 5, the values of cj,k, k = 0, 1, 2, ...,m; j = 1, 2, 3, 4 are found by solving a system of algebraic equations that is obtained depending on Tk(t) coefficients. 5. Numerical simulation To see how accurate the numerical scheme is, we run simulations on specific models in the range [0, 3] for the problem under study (1). The behavior of ψk(t) for k = 1, 2, 3, 4 is illustrated in Figures 1-5 at various parameter values p, m, ϱ. We examine the current model (1) with the subsequent parameter values [34]: θ1 = 0.2, θ2 = 0.02, λ1 = 0.3, λ2 = 0.4, d = 0.5, α1 = 6× 10−8, α2 = 3× 10−7, α3 = 1×10−7, α4 = 0.97, ρ = 0.01, κ = 2, β1 = 0.15, β2 = 0.7, β3 = 0.1, p = 0.5, σ = 0.1, p1 = 0.1, p2 = 0.2, p3 = 0.3, ϱ = 0.8. The I.Cs are ψi(0) = 0.2 for i = 1, 2, 3, 4. The numerical simulation of the examined model utilizing the specified methodology is illustrated in Figures 1-5. 1. Figure 1 presents the numerical solution for several quantities of the approximation order m = 5, 10, 15. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 9 of 16 2. Figure 2 illustrates the numerical solution for various values of ϱ = 0.8, 1.4, 2.0, 2.6. 3. Figure 3 illustrates the influence of p = 0.25, 0.5, 0.75, 1.0 on the numerical solution. 4. Figure 4 shows a comparison between the proposed method’s solution and the nu- merical solution found by the RK4 method [35], where the initial conditions and parameters were kept the same. 5. Figure 5 illustrates the residual error function (REF) [36] of the derived approxima- tion solution. The conduct of the numerical solution is contingent upon m, ϱ, p, as seen in Figures 1-3. Figures 2 and 3 indicate that the solution’s behavior aligns with the inherent influence of the parameters ϱ and p. Figures 4 and 5 indicate that the proposed strategy has been effectively utilized to address the problem under investigation. As a result, we can confirm that the disease behaved as expected. This means that we obtained an accurate simulation of the system that the relevant authorities can use to treat this deadly cancer. Figure 1. The approximate solution via various values of m. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 10 of 16 Figure 2. The approximate solution via various values of ϱ. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 11 of 16 Figure 3. The approximate solution via various values of p. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 12 of 16 Figure 4. Comparison of the solution obtained by the proposed method and RK4M. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 13 of 16 Figure 5. The REF of the obtained approximate solution. 6. Conclusions In this research, CHPM is applied to obtain numerical solutions for the BC with different initial solutions, and some parameters. By comparing the approximate solutions and the RK4 method ’s solution of the model under study, we were able to conclude that the approximate solutions obtained by applying the given technique are in excellent agreement with the RK4 method’s solution. Also through the resulting numerical results, we found to a large extent how effective this approach is in solving the problem under study and highlights the validity and potential of the proposed technique. Finally, this view of analytical and numerical solutions of dynamic variables was due to their being present and effectively influencing various models and fields of applied mathematics. Finally, the present study may contribute to providing more robust physical explanations for future theoretical and computational studies on the same topic. Acknowledgements The researchers wish to extend their sincere gratitude to the Deanship of Scientific Research at the Islamic University of Madinah for the support provided to the Post- Publishing Program. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 14 of 16 References [1] M. Fathoni, G. Gunardi, F. A. Kusumo and S. H. Hutajulu. Mathematical model analysis of breast cancer stages with side effects on the heart in chemotherapy patients, In Proceedings of the AIP Conference Proceedings, Yogyakarta, Indonesia, 29 July–1 August 2019; AIP Publishing: Melville, NY, USA, 2019. [2] Breast Cancer. Available online: https://www.who.int/news-room/fact- sheets/detail/breast-cancer (accessed on 15 July 2023). [3] C. E. De Santis, F. Bray, J. Ferlay, J. Lortet-Tieulent, B. O. Anderson and A. Jemal. International variation in female breast cancer incidence and mortality rates. Cancer Epidemiol Biomarkers Prev., 6:1495–1506, 2015. [4] M. Idrees, A. S. Alnahdi and M. B. Jeelani. Mathematical modeling of Breast Can- cer based on the Caputo-Fabrizio fractal-fractional derivative. Fractal Fract., 7:1–14, 2023. [5] M. Adel, M. M. Khader, T. A. Assiri andW. Kaleel. Numerical simulation for COVID- 19 model using a multidomain spectral relaxation technique. Symmetry, 15:1–12, 2023. [6] M. M. Khader, N. H. Sweilam, A. M. S. Mahdy and N. K. A. Moniem. Numerical simulation for the fractional SIRC model and influenza A. Appl. Math. Inf. Sci., 8:1029–1036, 2014. [7] M. Abaid Ur Rehman, J. Ahmad, A. Hassan, J. Awrejcewicz, W. Pawlowski, H. Karamti and F. M. Alharbi. The dynamics of a fractional-order mathematical model of Cancer Tumor Disease. Symmetry, 14:1–15, 2022. [8] M. M. Al-Shomrani and M. A. Abdelkawy. Numerical simulation for a fractional-order differential system of a Glioblastoma Multiforme and Immune system. Adv. Differ. Equ., 2020:516, 2020. [9] I. M. Batiha, A. Bataihah, Abeer A. Al-Nana, S. Alshorm, I. H. Jebril and A. Zraiqat. A numerical scheme for dealing with fractional initial value problems. Int. J. Innov. Comput. Inform. Control, 19:763–774, 2023. [10] F. Okose, M. T. Senel and R. Habbireeh. Fractional-order mathematical modeling of cancer cells-cancer stem cells-immune system interaction with chemotherapy. Math. Model. Numer. Simul. Appl., 1:67–83, 2021. [11] A. R. Alharbi, M. B. Almatrafi and Kh. Lotfy, Constructions of solitary traveling wave solutions for Ito integro-differential equation arising in plasma physics. Results in Physics, 19:103533, 2020. [12] A. M. S. Mahdy, Kh. Lotfy and A. A. El-Bary. Use of optimal control in studying the dynamical behaviors of fractional financial awareness model. Soft Computing, 26(7):3401–3409, 2022. [13] A. M. S. Mahdy, M. S. Mohamed, Kh. Lotfy, M. Alhazmi, A. A. El-Bary and M. H. Raddadi, Numerical solution and dynamical behaviors for solving fractional nonlinear Rubella ailment disease model. Results in Physics, 24(2):104091, 2021. [14] T. A. J. Al-Griffi and A. S. J. Al-Saif. Yang transform-homotopy perturbation method for solving a non-Newtonian viscoelastic fluid flow on the turbine disk. Z. Angew. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 15 of 16 Math. Mech., e202100116, 2022. [15] T. A. J. Al-Griffi and A. S. J. Al-Saif. Akbari-Ganji homotopy perturbation method for analyzing the pulsatile blood flow in tapered stenosis arteries under the effect of magnetic field together with the impact of mass and heat transfer. J. Comput. Appl. Mech., 53(4):543–570, 2022. [16] H. Aminikhah and M. Hemmatnezhad. An efficient method for quadratic Riccati differential equation. Commun. Nonlinear Sci. Numer. Simulat., 15:835–839, 2010. [17] R. Kumar, A. K. Singh and S. S. Yadav. New homotopy perturbation method for analytical solution of telegraph equation. Turk. J. Comput. Math. Educ., 12:2144– 2155, 2021. [18] K. Pal, V. G. Gupta, H. Singh and V. Pawar. Enlightenment of heat diffusion using new homotopy perturbation method. J. Appl. Sci. Eng., 27(3):2213–2216, 2023. [19] J. Wu, Y. Zhang, L. Chen and Z. Luo. A Chebyshev interval method for nonlinear dynamic systems under uncertainty. Math. Modell., 37:4578–4591, 2013. [20] Y. M. Hamada. A new accurate numerical method based on shifted Chebyshev se- ries for nuclear reactor dynamical systems. Sci. Technol. Nucl. Install, 15:ID7105245, 2018. [21] K. K. Ali, M. A. Abd El Salam, E. M. Mohamed, B. Samet, S. Kumar and M. S. Osman. Numerical solution for generalized nonlinear fractional integro-differential equations with linear functional arguments using Chebyshev series. Adv. Difference Equations, 1:1–23 2020. [22] F. Wang, Q. Zhao, Z. Chen and C. M. Fan. Localized Chebyshev collocation method for solving elliptic partial differential equations in arbitrary 2D domains. Appl. Math. Comput., 397:1–13, 2021. [23] Izadi, M. and Yüzbaşı, Ş. and Baleanu, D.. Taylor-Chebyshev approximation tech- nique to solve the 1D and 2D nonlinear Burger’s equations. Math. Sci., 16:459–471, 2022. [24] M. M. Khader and M. Adel. Modeling and numerical simulation for covering the fractional COVID-19 model using spectral collocation-optimization algorithms. Frac- tal Fract., 6:1–19, 2022. [25] A. d’Ornafrio, U. Ledzewicz H. Maurer and H. Schaettler. On optimal delivery of combination therapy for tumors. Math. Biosci., 222:13–26, 2009. [26] W. O. Kermarck and A. G. M. Kendrick. Contributions to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. A, 115:700-721, 1927. [27] C. Mufudza, S. Walter and E. T. Chiyaka. Assessing the effects of estrogen on the dynamics of chemo-virotherapy cancer. Comput. Math. Methods Med., 473572:1-15, 2012. [28] Z. Sabir, M. Munawar, M. A. Abdelkawy, M. A. Z. Raja, C. Ünlü, M. B. Jeelani and A. S. Alnahdi. Numerical investigations of the fractional-order mathematical model underlying Immune-Chemotherapeutic treatment for Breast Cancer using the neural networks. Fractal Fract., 6:1–16, 2022. [29] A. Yousef, F. Bozkurt and T. Abdeljawad. Mathematical modeling of the immune- chemotherapeutic treatment of breast cancer under some control parameters. Adva. M. Adel, M. M. Khader, M. Messaoudi / Eur. J. Pure Appl. Math, 18 (2) (2025), 6166 16 of 16 Differ. Equs., 696:1–25, 2020. [30] S. I. Oke, M. B. Matadi and S. S. Xulu. Optimal control analysis of a mathematical model for breast cancer. Math. Comput. Appl., 21:1–16, 2018. [31] J. C. Mason and D. C. Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, 2002. [32] M. H. Mudde, Chebyshev Approximation, University of Groningen, Netherlands, Fac- ulty of Science and Engineering, 2017. [33] H. C. Thacher. Conversion of a power to a series of Chebyshev polynomials. Commun. ACM, 7:181–182, 1964. [34] M. A. Hammad, I. H. Jebril, S. Alshorm, I. M. Batiha and N. A. Hammad. Numeri- cal solution for a fractional-order mathematical model of immune-chemotherapeutic treatment for breast cancer using the modified fractional formula. Int. J. Anal. Appl., 21:1–9, 2023. [35] M. M. Khader and Ali H. Tedjani. Using the modified decomposition method associ- ated with Mohand transforms for a numerical simulation of the mathematical model for Breast Cancer. European J. of Pure and Applied Mathematics, 18:1–13, 2025. [36] K. Parand and M. Delkhosh. Operational matrices to solve nonlinear Volterra- Fredholm integro-differential equations of multi-arbitrary order.Gazi University Jour- nal of Science, 29:895–907, 2016.