Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 427 https://internationalpubls.com Homotopy Analysis Method for the Approximate Solution of the SIRC Epidemic Model S. Geethamalini 1, S. Sangeetha2 , P. Venkataraman 3∗ 1School of Science and Humanities, Department of Mathematics, Sathyabama Institute of Science and Technology, Chennai-600119, Tamil Nadu, India. 2Department of Mathematics, SRM Institute of Science and Technology, Ramapuram-600089, Tamil Nadu,, India. 3Department of Medical Research, Faculty of Medicine and Health Sciences, SRM Institute of Science and Technology, Kattankulathur-603203, Tamil Nadu,, India. 1e-mails: 1geetha_malini@hotmail.com 2sangeets12@srmist.edu.in 3∗venkatap@srmist.edu.in Article History: Received: 22-10-2024 Revised: 06-12-2024 Accepted: 13-12-2024 Abstract: This paper investigates the Homotopy Analysis Method (HAM) as a means of approximating a solution to the SIRC epidemic model. This method allows the solution of the governing differential equation to be found as an infinite series with easily calculated components. The HAM uses an auxiliary parameter and a straightforward way to regulate and alter the region where the infinite series solution converges. The outcomes are displayed, and with just six terms, an extremely precise approximation solution can be achieved. Math. Subject Classification: 34G20; 34A34 Keywords:HAM, SIRC epidemic model, h-curve, Numerical Analysis 1. Introduction In the study and control of infectious diseases, mathematical models have grown to be crucial tools. In Kermack and McKendrick [1], one of the earliest models in epidemiology was presented to predict how a disease will spread. The total population is separated into three classes in this model: susceptible, infectious, and recovered, with the assumption being that it will remain stable over time. These models are hence known as SIR models. Cross-immune individuals (C) in the population have just recently been introduced in [2], they exist in a state that is between totally protected and unprotected (R). Because of this, the derived SIRC model considers transient partial immunity and might effectively characterise, say, influenza A. This study introduces and develops the homotopy analysis approach for approximately solving the SIRC model. Parameters and variables are presented in Table 1, and the model is in Table 1, and the model is (1 )S S SI C   • = − − + ( )I SI CI I    • = + − + (1 ) ( )R CI I R     • = − + − + (1) ( )C CI C R    • = − − + + Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 428 https://internationalpubls.com Table I. Values of parameters and variables Parameters and Variables Meaning  Death amount in every division presumed to be equal to the amount of new born in the population  Amount of re-susceptibility of the cross-immune population  Rate of contact  Average probability of reinfection of cross-immune individuals  Regaining amount of the infected population  Amount at which the regaining population to the cross-immune population and from fully immunized to partial immunity S Susceptible I Infected R Recovered C Cross-immune Rihan et al., [3] studied the fractional SIRC model with salmonella bacterial infection. Amjad et al., [4] have studied the numerical simulation of SIRC model of fractional order derivative. Geethamalini et al. [5] have published numerical and analytical study of sirc epidemical model using HPM.Ghoreishi et., al. [6] have studied the HAM for solving for CD4+ T-cells of HIV infection. The authors was published a semi-analytical solutions of mathematical models in EIAV (Equine Infectious Anemia Virus) infection using HAM. Also the authors in [7–9] were established an approximate solutions and dynamical analysis of an EIAV infection . Alijhani et al. [10] were studied the numerical solution of Fractional-order HIV model using HAM. Naik et al. [11–13] have studied the utilizing the homotopy analysis method to estimate the approximate analytical solution of the HIV viral dynamic model of CD4+ T-cells. (Odibat and Baleanu, [14] studied the linearization based approach of HAM for nonlinear time fractional parabolic PDEs. Saad et al. [15] studied the exact solutions for time fractional Burger’s equations using HAM. Yepez and Gomez, [16] was introduced, an updated explanation of the Caputo Fabrizio fractional-order derivative, including with its applications to the multistep HAM Deniz, [17 ]have studied the Semi-analytical approach for solving a model for hiv infection. In recent years, different nonlinear systems of differential equations inmathematics and the sciences have been solved using HAM as a solution technique [18–21]. The HAM was first developed by [22–24]. Chioma et al. [25] studied the Application of Homotopy Analysis method for solving an SEIRS epidemic Model. Abdi-Rahim et al. [26] studied the Analytical Study of Fractional Epidemic Model via Natural Transform Homotopy Analysis Method. The smoking pandemic model with fractional order was investigated by Veeresha et al. [27] using a modified homotopy analysis transforms approach. Bakare et al. [28] studied the Interval-based uncertain SIR epidemic model numerically solved using HAM. Duarte et al. [29] studied the Chaos analysis of SIR epidemic Method. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 429 https://internationalpubls.com 2. Solution of the SIRC model by the HAM In order to develop the HAM solutions to (1), first we select 0 0 0 0(0) , (0) , (0) , (0) . (2)S S I I R R C C= = = = The auxiliary linear operators 1 2,,L L , and 3L 4L are selected as 1 ( , ) [ ( , )] , dS t s L S t s dt = 2 ( , ) [ ( , )] , dI t s L I t s dt = 3 ( , ) [ ( , )] . dR t s L R t s dt = 4 ( , ) [ ( , )] . dC t s L C t s dt = Which satisfy the following properties: ( ) 0,i iL C = where ( 1, 2,3, 4)iC i = are integral constants. Describe 1[ , , , ] (1 ) ,N S I R C S S SI C   • = − − + − 2[ , , , ] ( ) ,N S I R C I SI CI I    • = − − + + 3[ , , , ] (1 ) ( )N S I R C R CI I R     • = − − − + + 4[ , , , ] ( )N S I R C C CI C R    • = + + + − Introduce nonzero auxiliary parameter h and nonzero auxiliary function using Liao's definitions. The embedding parameter is H (t) and s ∈ [0, 1]. We create the zero-order deformation equations using this. 1 0 1 1 1(1 ) [ ( ; ) ( )] ( ) [ , , , ], (3)s L S t s S t sh H t N S I R C− − = 2 0 2 2 2(1 ) [ ( ; ) ( )] ( ) [ , , , ], (4)s L I t s I t sh H t N S I R C− − = 3 0 3 3 3(1 ) [ ( ; ) ( )] ( ) [ , , , ], (5)s L R t s R t sh H t N S I R C− − = 4 0 4 4 4(1 ) [ ( ; ) ( )] ( ) [ , , , ], (6)s L C t s C t sh H t N S I R C− − = clearly, when s=0 and s=1, then 0( ;0) ( ), ( ;1) ( ),S t S t S t S t= = 0( ;0) ( ), ( ;1) ( ),I t I t I t I t= = Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 430 https://internationalpubls.com 0( ;0) ( ), ( ;1) ( ).R t R t R t R t= = 0( ;0) ( ), ( ;1) ( ).C t C t C t C t= = When a result, as s the embedding parameter increases from 0 to 1, the solutions S (t; s), I (t; s) ,R (t; s) & C(t;s) varies continuously from S0 (t), I0 (t) , R0 (t) & C0((t) to the exact solution S (t), I (t), R(t) & C (t). Using Taylor’s series 0 1 ( ; ) ( ) ( ) , (7)i i i S t s S t S t s  = = + 0 1 ( ; ) ( ) ( ) , (8)i i i I t s I t I t s  = = + 0 1 ( ; ) ( ) ( ) , (9)i i i R t s R t R t s  = = + 0 1 ( ; ) ( ) ( ) , (10)i i i C t s C t C t s  = = + where 0 1 ( ; ) | , ! i i si S t s S i s =  =  0 1 ( ; ) | , ! i i ss I t s I i s =  =  0 1 ( ; ) | , ! i i si R t s R i s =  =  0 1 ( ; ) | , ! i i si C t s C i s =  =  If h1, h2, h3,h4, H1(t),H2(t), H3(t) and H4(t) are selected, and the series converges at p=1. 0 1 ( ) ( ) ( ),i i S t S t S t  = = + 0 1 ( ) ( ) ( ),i i I t I t I t  = = + 0 1 ( ) ( ) ( ).i i R t R t R t  = = + 0 1 ( ) ( ) ( ).i i C t C t C t  = = + The ith-order deformation equations are obtained by differentiating (3)–(6) ‘i’ times with regard to s, allocating by i! and setting s = 0. 1 1 1, 1[ ( ) ( )] ( ( )),i i i i iL S t S t hR S t − −− = (11) 2 1 2, 1[ ( ) ( )] ( ( )),i i i i iL I t I t hR I t − −− = (12) 3 1 3, 1[ ( ) ( )] ( ( )),i i i i iL R t R t hR R t − −− = (13) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 431 https://internationalpubls.com 4 1 4, 1[ ( ) ( )] ( ( )),i i i i iL C t C t hR C t − −− = (14) where 1 1 1, 1 1 1 0 ( ) ( ) ( ) ( ) ( ) ( ) (1 ) , i i i j i j i i i j dS t R t S t I t S t C t dt      − − − − − − = = + + − − − 1 1 1 2, 1 1 1 0 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ), i i i i j i j j i j i j j dI t R t S t I t C t I t I t dt     − − − − − − − − = = = − − + +  1 1 3, 1 1 1 0 ( ) ( ) (1 ) ( ) ( ) ( ) ( ) ( ), i i i j i j i i j dR t R t C t I t I t R t dt      − − − − − − = = − − − + + 1 1 4, 1 1 1 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ), i i i j i j i i j dC t R t C t I t C t R t dt     − − − − − − = = + + + − & 1 1 . 0 1 i i i    =     When i is higher than or equal to 1, ith-order deformation (11)–(14) becomes 1 1, 0 ( ) ( ) ( ) , t i i i iS t S t h R d  −= +  1 2, 0 ( ) ( ) ( ) , t i i i iI t I t h R d  −= +  1 3, 0 ( ) ( ) ( ) . t i i i iR t R t h R d  −= +  1 4, 0 ( ) ( ) ( ) . t i i i iC t C t h R d  −= +  3. Numerical Simulations Take into account the values below for the numerical outcomes [3]. 0 0 0 00.3, 0.5, 0, 0.6S I R C= = = = 1.3, 0.09, 0.1, 0.05, 0.36, 0.9     = = = = = = We obtain the sixth-order for using the Mathematica software. ( ), ( ), ( ) ( )S t I t R t and C t were obtained, and are given below 2 3 4 5 6 2 2 3 2 4 2 5 2 6 2 3 3 4 3 5 3 6 3 4 4 5 4 ( ) 0.3 0.612 1.53 2.04 1.53 0.612 0.102 0.550575 1.4682 1.65173 0.88092 0.183525 0.667932 1.50285 1.20228 0.333966 0.0627671 0.100427 0.0 S t ht h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t = + + + + + + − − − − − − − − − + + + 6 4 5 5 6 5 6 6 418447 0.0461811 0.0384842 0.00050418 ... (15) h t h t h t h t + + − + Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 432 https://internationalpubls.com 2 3 4 5 6 2 2 3 2 4 2 5 2 6 2 3 3 4 3 5 3 6 3 4 4 5 4 6 ( ) 0.5 1.926 4.815 6.42 4.815 1.926 0.321 1.03131 2.75016 3.09393 1.6501 0.34377 1.63923 3.68826 2.95061 0.819613 0.2508 0.40128 0.1672 I t ht h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t h = − − − − − − − − − − − + − − + + + + 4 5 5 6 5 6 60.115992 0.09666 0.00467747 ... (16) t h t h t h t − + − + 2 3 4 5 6 2 2 3 2 4 2 5 2 6 2 3 3 4 3 5 3 6 3 4 4 ( ) 0.001 1.314 3.285 4.38 3.825 1.314 0.219 0.511335 1.36356 1.53401 0.81813 0.170445 0.344179 0.774402 0.619522 0.172089 0.0790303 0.126448 R t ht h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t h = − − − − − − + + + + + + + + + − − 5 4 6 4 5 5 6 5 6 6 0.0526868 0.0173387 0.0144489 0.000654392 ... (17) t h t h t h t h t − − − + + 2 3 4 5 6 2 2 3 2 4 2 5 2 6 2 3 3 4 3 5 3 6 3 4 4 5 4 6 4 ( ) 0.6 2.844 7.11 9.48 7.11 2.844 0.474 1.09484 2.9196 3.28455 1.7516 0.36495 1.3145 2.95763 2.3661 0.65725 0.23452 0.375233 0.156347 C t ht h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t h t = + + + + + + + + + + + − − − − − − − 5 5 6 5 6 60.0871498 0.0726248 0.00452726 ... (18)h t h t h t + + + + 4. Discussion We found the solution from (11)-(14) which contain "h" demonstrate a simple method Liao suggests for controlling and adjusting curves to validate series solutions converge. The graphs of the 5th and 6th term approximations of S, I, R &C is shown in Figures 1–5. From these curves it shows the the horizontal axis forms a valid region of ’”h” Table 2 contains a list of the suitable regions. To get the ideal values for h, an error analysis is performed. We enter Eqs (15) through (18) into (1) and obtain the corresponding residual functions: 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 ( ; ) ( , , , ; ) (1 ) ( ; ) ( ; ) ( ; ) (19) ( ; ) ( , , , ; ) ( ; ) ( ; ) ( ; ) ( ; ) ( ) ( ; ) (20) ( ; ) ( , , , ; ) (1 ) ( ; ) ( ; ) S S I C I S I C I I R C I d t h ER S I R C h S t h t h t h dt d t h ER S I R C h t h t h t h t h t h dt d t h ER S I R C h t h t h dt                   = − − + − = − − + + = − − − 3 3 4 4 4 4 4 4 4 ( ; ) ( ) ( ; ). (21) ( ; ) ( , , , ; ) ( ; ) ( ; ) ( ) ( ; ) ( ; ). (22) I R C C I C R t h t h d t h ER S I R C h t h t h t h t h dt            + + = + − + + Based on the following [27-29], for the 6th order approximation, we estimate the square residual error to be Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 433 https://internationalpubls.com 1 2 1 1 1 0 1 2 2 2 2 0 1 2 3 3 3 0 1 2 4 4 4 0 ( ) ( ( , , , ; )) , (23) ( ) ( ( , , , ; )) , (24) ( ) ( ( , , , ; )) , (25) ( ) ( ( , , , ; )) . (26) RS h ER S I R C h dt RI h ER S I R C h dt RR h ER S I R C h dt RC h ER S I R C h dt = = = =     Table 2. Figures (1)–(8) shows the values of .h ( ) 1.4 0.4S t h−   − ( ) 1.3 0.6I t h−   − ( ) 1.3 0.6R t h−   − ( ) 1.4 0.5C t h−   − Figure 1: The h-curves of S′(0) obtained by the fifth-order and sixth-order approximations of the HAM Figure 2: The h-curves of I′(0) obtained by the fifth-order and sixth-order approximations of the HAM. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 434 https://internationalpubls.com Figure 3: The h-curves of R′(0) obtained by the fifth-order and sixth-order approximations of the HAM. Figure 4: The h-curves of C′(0) obtained by the fifth-order and sixth-order approximations of the HAM. (a) (b) (c) (d) Figure 5: (a) Various h = −0.9, h = −1, h = ideal, h = −1.1 and t ∈ (0, 1) the residual errors function of Eq.(19). (b) Various h = −0.9, h = −1, h = ideal, h = −1.1 and t ∈ (0, 1) the residual errors function of Eq.(20). © Various h = −0.9, h = −1, h = ideal, h = −1.1 and t ∈ (0, 1), the residualerrors function Eq. of (21). (d) Various h = −0.9, h = −1, h = ideal, h = −1.1 and t ∈ (0, 1), the residual errors function Eq. of (22) Values of h1, h2, h3 and h4 for which ￣RS(h1), ￣RI(h2), ￣RR(h3) & ￣RC(h4) are minimum. Then Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 435 https://internationalpubls.com ** * * 31 2 4 1 2 3 3 ( )( ) ( ) ( ) 0, 0, 0, 0.. dRR hdRS h dRI h dRC h dh dh dh dh = = = = The optimal values h1,h2 , h3 & h4 for all of the cases considered are obtained as * * * * 1 2 3 40.856097, 1.0557, 0.913549, 1.04116.h h h h= − = − = − = − Table 3 shows the minimal values of 1 2 3 4RS (h ), RI (h ),RR (h ) & RC (h ) for the optimal values of h1, h2 h3 and h4 . We estimated errors for various t in (0, 1) which is listed in Table 4. This demonstrates that the HAM provides us with a close approximation to the solution for the SIRC model (1). The residual errors of t in (0,1) and different h are shown in Figures 5(a,b,c &d). Figure 6: The optimum and Minimum value of S(t). Figure 7 : The optimum and Minimum value of I(t). Figure 8: The optimum and Minimum value of R(t). Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 436 https://internationalpubls.com Figure 9: The optimum and Minimum value of C(t). Table 3: The minimum values of 1 2 3 4RS (h ), RI (h ),RR (h ) & RC (h )     (see Figures 6-9) h* Minimum value 8 1( ) 0.856097 8.38046 10RS h −−  8 2( ) 1.0557 9.75647 10RI h −−  9 3( ) 0.913549 2.79827 10RR h −−  9 4( ) 1.04116 9.28625 10RC h −−  Table 4: The residual errors 1 2 3 4ER , ER , ER & ER for various (0,1).t t * 1 1( , , , ; )ER S I R C h * 2 2( , , , ; )ER S I R C h * 3 3( , , , ; )ER S I R C h * 4 4( , , , ; )ER S I R C h 0.0 79.05773 10− 99.58597 10− 89.14248 10− 92.30479 10− 0.1 61.52525 10− 95.94778 10− 71.69604 10− 88.5875 10− 0.2 67.34594 10− 66.4794 10− 75.70379 10− 63.30688 10− 0.3 63.06653 10− 53.89008 10− 62.02822 10− 51.80537 10− 0.4 55.53734 10− 41.19367 10− 51.20506 10− 55.20511 10− 0.5 41.66688 10− 42.54236 10− 53.24383 10− 41.03328 10− 0.6 43.26493 10− 44.1493 10− 56.11003 10− 41.53968 10− 0.7 44.74678 10− 45.24149 10− 58.70219 10− 41.71343 10− 0.8 44.78428 10− 44.52458 10− 58.63747 10− 41.20174 10− 0.9 41.09547 10− 52.87077 10− 51.88407 10− 41.26712 10− 1 49.76499 10− 49.33951 10− 41.75798 10− 41.80173 10− 5. Conclusion In order to solve an epidemic SIRC model, the homotopy analysis method has been effectively developed and utilised in this study. The auxiliary parameter h, which is present in the HAM solution, provides an easy method for adjusting and controlling the convergence region of the resulting infinite series. The outcomes demonstrate that the HAM is an exact and effective method for determining the approximation. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 437 https://internationalpubls.com Acknowledgements: The authors are grateful to anonymous referees for their excellent suggestions, which greatly improved the presentation of the paper. Declarations Competing interests The authors declare that there is no conflict of interests. Refrences [1] Kermack, W.O., McKendrick, A.G.: Contributions to the mathematical theory of epidemics: II. The problem of endemicity. Bull Math Biol. 53(1-2), 57-87, (1991). [2] Casagrandi, R., Bolzoni, L., Levin, et al.: The SIRC model and influenza A, Math Biosci. 200(2), 152-169, (2006). [3] Rihan, F.A., Alsakaji, H.J., Rajivgandhi, C.: Stochastic SIRC epidemic model with timedelay for COVID-19. Adv Differ Equ. 2020(1), 502, (2020). [4] Amjad, A., Muhammaad, Y.K., Gouhar, A.: Qualitative theory and numerical simulation of SIRC model corrresponding to nonlocal fractional order derivative. Comm Nonlinear Anal. 2, 75-91, (2020). [5] Geethamalini,S., Sangeetha,S., Anandhababu,D., Venkataraman,P.: A Numer ical and Analytical Study of Sirc Epidemic Model Using HPM. J. of Com put.Ana.and Appli. 33(7), 1097-1104(2024). [6] Ghoreishi, M., M Ismail, A.I.B., Alomari, A.K.: Application of the homotopy analysis method for solving a model for HIV infection of CD4+ T-cells, Math Comput Model.54(11-12), 3007-3015, (2011). [7] Geethamalini, S., Balamuralitharan, S.: Semi-analytical solutions by homotopy analysis method for EIAV infection with stability analysis. Adv Differ Equ. 356, 1-14, (2018). [8] Balamuralitharan, S., Geethamalini, S.: Solutions of the epidemic of EIAV in fection by HPM. J Phys Conf Ser. 1000(1), 1-7, (2018). [9] Geethamalini, S., Balamuralitharan, S.: Dynamical analysis of EIAV infection with cytotoxic T lymphocyte immune response delay. RINAM. 2 100025, (2019). [10] Aljhani, S., Noorani, M.S., Alomari, A.K.:Numerical Solution of Fractional Order HIV Model Using Homotopy Method. Discrete Dyn Nat Soc. 2020, 1-13 (2020). [11] Naik, P.A., Jian, Zu., Ghoreishi, M.: Estimating the approximate analytical solution of HIV Viral dynamic model by using homotopy analysis method. Chaos Solitons Fractals.131,109500, (2019). [12] Naik, P.A., Ghoreishi, M., Jian, Zu.: Approximate Solution of a nonlinear fractional- order HIV model using homotopy analysis method. Int. J. Numer. Anal. Model.19(1) 52-84, (2022). [13] Naik, P.A., Ghoreishi, M.: Stability analysis and approximate solution of SIR epidemic model with Crowley-Martin type functional response and Holling type II treatment rate by using homotopy analysis method. J. Appl. Anal. Comput. 10(4) 1482-1515 (2020). [14] Odibat, Z., Baleanu, D. A linearization-based approach of homotopy analysis method for non-linear time-fractional parabolic PDEs. Mathematical Methods in Applied Sciences, 42(18) 7222–7232,(2019). [15] Saad, K.M., AL-Shareef, E H.F., Alomari, A. K., Baleanu, D.: On exact solutions for time-fractional Korteweg-de Vries and Korteweg-de Vries-Burger’s equa tions using homotopy analysis transform method. Chinese Journal of Physics, 63, 149–162,(2020). [16] Yepez-Martnez, H., Gomez-Aguilar, J. F. A new modified definition of Caputo Fabrizio fractional- order derivative and their applications to the multistep ho motopy analysis method. Journal of Computational Applied and Mathematics, 346, 247–260,(2019). [17] Deniz S.: Semi-analytical approach for solving a model for HIV infection of CD4+ T-cells. TWMS J. of Apl. and Eng. Math. 11(1),273-281 (2021). [18] Awawdeh, F., Adawi, A., Mustafa, Z. : Solutions of the SIR models of epidemics using HAM. Chaos Solitons Fractals. 42, 3047-3052, (2009). [19] Momoh, A.A., Ibrahim, M.O., Tahir, A., Adamu, I.I.: Application of Homotopy Analysis Method for Solving SEIR models of Epidemics. Nonlinear Differ. Equ. Appl. 3(2), 53-68, (2015). [20] Nirmala, P., Subramanian, S. P.: SEIR Model of Seasonal Epidemic Diseases using HAM. Appl Appl Math. 10(2), 1066-1081, (2015). Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 438 https://internationalpubls.com [21] Khan, H., Mohapatra, R.N., Vajravelu, K., Liao, S.: The explicit series solution of SIR and SIS epidemic models. Appl Math Comput. 215(2), 653-669, (2009). [22] Liao, S. J.: Beyond perturbation: introduction to the homotopy analysis method, CRC Press, Chapman and Hall, Boca Raton (2003). [23] Liao, S. J.: Comparison between the homotopy analysis method and homotopy perturbation method. Appl Math Comput. 169(2),1186-1194, (2005). [24] Liao, S. J.: An optiomal homotopy-analysis approach for strongly nonlinear differential equation. Commun Nonlinear Sci Numer Simul. 15(8), 2003-2016, (2010). [25] Chioma, I.S., Ugonna, E.G., Michael, U.O., et.al.: Application of Homotopy Analysis method for solving an SEIRS epidemic Model. Math.Model.Appl. 4(3),36-48, (2019). [26] Abdi-Rahim,H.R., Zayed, M., Ismail, G.M. . Analytical Study of Fractional Epi demic Model via Natural Transform Homotopy Analysis Method. Symmetry, 14(8), 1-18,(2022). [27] Veeresha, P., Prakasha, D.G., Baskonus, H.M.: Solving smoking epidemic model of fractional order using a modified homotopy analysis transform method.Matheatical Sciences, 13, 115–128,(2019). [28] Bakare, E.A., Chakraverty, S., Potucek, R.: Numerrical solution of an interval based uncertain SIR epidemic model by Homotopy Analysis Method. Axioms. 10(2), 1-19, (2021). [29] Duarte, J., Januario, C., Martins, N., etal.: Chaos analysis and explicit series solutions to the seasonally forced SIR epidemic model. J. Math. Biol. 78, 2235- 2258, (2019).