1 Introduction

One of the most important diseases faced by human beings ever is tuberculosis (TB). TB is a spreadable, airborne bacterial infection which is caused by mycobacterium tuberculosis. This bacterium normally affects the lungs (pulmonary tuberculosis). This bacterium may also affect several other systems like kidneys, the brain, the lymphatic system, the central nervous system, spinal cord, etc. The presence of TB disease has been found since ancient times in various civilizations like Egypt, China, Roman, etc. (see [1]). One-third of the world population at the present time is infected with TB, and the number of infectious individuals increases at a rate of one per second [2]. The aforementioned disease was among the top ten death causes around the globe in the year 2015, where about 10.4 million individuals were infected from it. In the same year, about 1.8 million infectious individuals lost their lives from various diseases including 0.4 million with HIV and TB. 60% of the tuberculosis cases around the globe were concentrated in the six countries (Pakistan, India, Nigeria, China, Indonesia, and South Africa) [3]. Dye [4] gave some information that the major cause of death worldwide, in particularly in Sub-Saharan Africa, is due to TB and HIV. Further, HIV epidemic is a serious threat in many countries of the world. It is a clear evidence that worldwide children are protected from the early infection of the disease by vaccination like Bacillus Calmatte–Guerine(BCG) [5]. Therefore, detection and treatment of latent TB by modern therapy have been used recently to prevent the breakdown of rate of spread of the disease as only the members of the infectious class can transmit the disease to others.

Worldwide various procedures and methods have been used to understand the cause and control of these diseases in society. One of the powerful tools is the mathematical modeling through which we understand the dynamics of various disease transmission and suggest procedure how to control them in society. The mentioned area originated during 1927 for the first time. A variety of models have been developed and studied so far (refer to [610]). In this regard a five-compartment model for TB has also been constructed in [11] as follows:

$$ \textstyle\begin{cases} \dot{\mathscr{M}}(t)=\theta \rho -(\alpha +\mu ) \mathscr{M}, \\ \dot{\mathscr{S}}(t)=(1-\theta )\rho +\alpha \mathscr{M}-\beta { \mathscr{SI}}-\rho {\mathscr{S}}, \\ \dot{\mathscr{L}}(t)=\beta \mathscr{SI}-(\sigma +\tau +\mu ) \mathscr{L}, \\ \dot{\mathscr{I}}(t)=\tau \mathscr{L}- (\gamma +\mu +\delta )\mathscr{I}, \\ \dot{\mathscr{R}}(t)=\sigma \mathscr{L}+\gamma \mathscr{I}-\mu \mathscr{R}. \end{cases} $$
(1)

In the above model the whole populace is categorized into five classes: the immunized class \(\mathscr{M}\), the susceptible class \(\mathscr{S}\), the infected latently class \(\mathscr{L}\), the infectious class \(\mathscr{I}\), and the recovered class \(\mathscr{R}\). Parameters of the model under consideration are explained as follows: The constant of recruitment is represented by the symbol ρ, θ denotes the immunized portion at birth, α represents the rate of weaning off the vaccine, the natural death rate is denoted by the symbol μ, β represents the tuberculosis contraction rate, the successful treatment of infectious latent is denoted by the letter σ, the symbol τ is the rate of breakdown of latent TB into infectious TB, the successful cure of infectious TB individuals and the death resulting from the disease are respectively denoted by the symbols γ, δ.

Usually an integer order derivatives do not explore the dynamics of real world problems related to biology and physics well. In order to overcome this deficiency, fractional calculus has been given attention for the last few decades. Also we know that fractional calculus is increasingly used by mathematicians for mathematical modeling. Derivatives and integrals of noninteger order may be defined by a number of ways. Some well-known definitions are those given by Riemann and Liouville [12], Caputo, etc. (see [13]). The mentioned derivatives involve kernel of singular type. Frequently these two definitions have been increasingly used since fractional differential operator is in fact a definite integral operator for whom the definition of kernel is not unique or not regular. Further due to high degree of freedom in derivative of arbitrary order, researchers have given much attention to studying applied problems under these concepts. In this sense very recently some authors replaced the singular kernel by some nonlocal nonsingular and produced new definitions. Hence in 2015 Caputo and Fabrizio replaced the singular kernel by exponential kernel in the usual Caputo derivative and called it Caputo–Fabrizio fractional derivative (abbreviated as CFFD); for details, see [1416]. Therefore various researchers investigated different problems of applied nature under this concept. In various cases the mentioned derivative has produced significant results as compared to other forms of derivatives (see [1719]). The Caputo–Fabrizio derivative omits singular kernel by exponential kernel and hence makes the concerned differential operator nonlocal. Conventional fractional derivatives contain singular kernel which sometimes causes difficulty in explanation of some characteristics of various materials. To overcome this, Caputo and Fabrizio introduced a new definition of fractional integral and derivative which involves exponential kernel instead of singular one. Various studies can be found in the literature that have focused on the Caputo–Fabrizio fractional order derivatives, see, for instance, [2024]. Further, in various papers related to thermal sciences, the mentioned derivative has been proved powerful as compared to other type, we refer to [2530]. Keeping these points in mind and the nonsingular nature of the proposed derivative, we investigate the considered model for existence and analytical results.

Now the question is how to treat the problems involving derivative of fractional orders. For this need researchers have successfully updated the usual tools and methods to handle differential equations of fractional order (FODEs). Usual perturbation techniques and decomposition methods were greatly utilized to deal with ordinary FODEs. Also, for the mentioned problems, Adomian decomposition coupled with some integral transforms has been used very well (see [3134]). On the other hand, since the FODEs involving new type derivatives are very rarely used, very frequently authors [35] established some algorithms to handle such FODEs containing CFFD.

Hence we investigate the model given in (1) under CFFD as follows:

$$ \textstyle\begin{cases} {} _{0}^{CF}D_{t}^{\eta } \mathscr{M}(t)=\theta \rho -(\alpha +\mu ) \mathscr{M}, \\ {} _{0}^{CF}D_{t}^{\eta } \mathscr{S}(t)=(1-\theta )\rho +\alpha \mathscr{M}-\beta \mathscr{SI}-\rho \mathscr{S}, \\ {} _{0}^{CF}D_{t}^{\eta } \mathscr{L}(t)=\beta \mathscr{SI}-(\sigma + \tau +\mu )\mathscr{L}, \\ {} _{0}^{CF}D_{t}^{\eta } \mathscr{I}(t)=\tau \mathscr{L}-(\gamma +\mu + \delta )\mathscr{I}, \\ {} _{0}^{CF}D_{t}^{\eta } \mathscr{R}(t)=\sigma \mathscr{L}+\gamma \mathscr{I}-\mu \mathscr{R}. \end{cases} $$
(2)

We study model (2) subject to the biologically feasible initial conditions

$$\begin{aligned}& \mathscr{M}(0)=\mathbf{N}_{1}\geq 0,\qquad \mathscr{S}(0)= \mathbf{N}_{2} \geq 0, \\& \mathscr{L}(0)=\mathbf{N}_{3} \geq 0, \\& \mathscr{I}(0)=\mathbf{N}_{4}\geq 0, \qquad \mathscr{R}(0)= \mathbf{N}_{5} \geq 0. \end{aligned}$$

Initially we establish some conditions about the existence of solutions for model (2) by using some fixed point results like Banach and Krasnoselskii. After that, by using the considered tool of “Laplace Adomian decomposition method (LADM)” for \(\eta \in (0,1]\), we compute semi-analytical results. Finally, the approximate results are presented via graphs.

2 Preliminaries

In this section of the manuscript, we present some fundamental definitions as follows.

Definition 1

([15])

Let \(\varphi \in \mathcal{H}^{1}(a,b)\), \(b>a\), \(\eta \in (0,1)\), then the CFFD is given as

$$ {}_{0}^{CF}D_{t}^{\eta } \varphi (t)=\frac{\mathcal{K(\eta )}}{1-\eta } \int _{\alpha }^{t} \varphi ^{\prime }(\theta ) \exp \biggl[- \frac{t-\theta }{1-\theta } \biggr]\,d\theta , $$
(3)

where \(\mathcal{K(\eta )}\) in (3) is the normalization function with \(\mathcal{K}(1)=\mathcal{K}(0)=1\). If the function failed to exist in \(\mathcal{H}^{1}(a,b)\), then the above derivative can be reformulated as

$$ {}_{0}^{CF}D_{t}^{\eta }\varphi (t)= \frac{\mathcal{K(\eta )}}{1-\eta } \int _{\alpha }^{t} \varphi (t)- \varphi (\theta ) \exp \biggl[- \frac{t-\theta }{1-\theta } \biggr]\,d\theta . $$

Definition 2

([36])

Let \(\eta \in (0,1]\), then the integral of fractional order η of the function φ is

$$\begin{aligned} {}_{0}^{CF}I_{t}^{\eta }\varphi (t)= \frac{(1-\eta )}{\mathcal{K}(\eta )} \varphi (t) +\frac{\eta }{\mathcal{K}(\eta )}\varphi (t) \int _{0}^{t} \varphi (\theta )\,d\theta . \end{aligned}$$

Definition 3

([35, 37])

The Laplace transform of CFFD \({}_{0}^{CF}D_{t}^{\eta }\) of \(M(t) \) is given as

$$ \mathscr{L} \bigl[{}_{0}^{CF}D_{t}^{\eta } M(t) \bigr]= \frac{s \mathscr{L}[M(t)]-M(0)}{s+\eta (1-s)}, \quad s\geq 0, \eta \in (0,1]. $$

2.1 Equilibrium points and the basic reproduction number

Before proceeding further, we consider it advantageous to find the equilibrium points and the basic reproduction number of the model under consideration. There are two types of possible equilibrium points of the model. The first one is the point where there is no disease in the community, i.e., the disease-free equilibrium point. This is obtained by setting the right-hand side of each equation in the model to zero along with \(\mathscr{L}=\mathscr{I}=\mathscr{R}=0\). Solving the system then gives \(\mathscr{M}^{0}=\frac{\theta \rho }{\alpha +\mu }\) and \(\mathscr{S}^{0}=\frac{\alpha +\mu -\mu \theta }{\alpha +\mu }\). Thus the disease-free equilibrium point of the model under investigation is given by \(\mathscr{E}^{0}= (\frac{\theta \rho }{\alpha +\mu }, \frac{\alpha +\mu -\mu \theta }{\alpha +\mu },0,0,0 )\).

To find the basic reproduction number, we consider only the infectious classes of the model. Let \(\mathbf{{\mathit{{V}}}= (\mathscr{L},\mathscr{I} )^{T}}\), by the help of the given model one can write d V d t =FV= [ β S I 0 ] [ ( σ + τ + μ ) L τ L + ( γ + μ + σ ) I ] . The Jacobian matrices of \(\mathscr{F}\) and \(\mathscr{V}\) are given by F= [ 0 β S 0 0 0 ] and V= [ σ + τ + μ 0 τ γ + μ + δ ] . The inverse matrix of \(\mathcal{V}\) is given by V 1 = [ 1 σ + τ + μ 0 τ 1 γ + μ + δ ] . Hence the next generation matrix \(\mathcal{F}\mathcal{V}^{-1}\) is calculated as

F V 1 = [ β S 0 τ ( γ + μ + δ ) ( σ + τ + δ ) β S 0 γ + μ + δ 0 0 ] .
(4)

The spectral radius of the next generation matrix (4) gives the threshold quantity \(R_{0}\) [38]. Thus

$$ R_{0}= \frac{\beta \tau (\alpha +\mu -\mu \theta )\tau }{(\alpha +\mu )((\gamma +\mu +\delta )(\sigma +\tau +\delta )}. $$

Three-dimensional plots of the basic reproduction number \(R_{0}\) versus different parameters in the model under consideration are depicted in Fig. 1. This quantity plays the key role in stability analysis and in finding conditions for the said purpose.

Figure 1
figure 1

Three-dimensional plot of the basic reproduction number versus different parameters involved in the model under consideration. The parametric values are given in Table 1 at the end

3 Existence and uniqueness results for tuberculosis disease model of fractional order

In the following we derive existence results related to our model (2) exploiting the so-called fixed point theorem due to Banach. To proceed further, we first of all define the functions given below:

$$ \textstyle\begin{cases} f_{1}(t,M,S,L,I,R)=\theta \rho -(\alpha +\mu )M, \\ f_{2}(t,M,S,L,I,R)=(1-\theta )+\alpha {M}-\beta {SI}-\rho {S}, \\ f_{3}(t,M,S,L,I,R)=\beta {SI}-(\sigma +\tau +\mu ){L}, \\ f_{4}(t,M,S,L,I,R)=\tau {L}-(\gamma +\mu +\delta ){I}, \\ f_{5}(t,M,S,L,I,R)=\sigma {L}+\gamma {I}-\mu {R}, \end{cases} $$
(5)

where

$$ M(0)=\mathbf{N}_{1},\qquad S(0)=\mathbf{N}_{2},\qquad L(0)= \mathbf{N}_{3},\qquad I(0)= \mathbf{N}_{4},\qquad R(0)= \mathbf{N}_{5}. $$

So our problem becomes

$$ \textstyle\begin{cases} {} _{0}^{CF}D_{t}^{\eta }M(t)= f_{1}(t,M,S,L,I,R)=\theta \rho -(\alpha + \mu )M, \\ {} _{0}^{CF}D_{t}^{\eta } S(t)= f_{2}(t,M,S,L,I,R)=(1-\theta )+\alpha {M}- \beta {SI}-\rho {S}, \\ {} _{0}^{CF}D_{t}^{\eta } L(t)= f_{3}(t,M,S,L,I,R)=\beta {SI}-(\sigma + \tau +\mu ){L}, \\ {} _{0}^{CF}D_{t}^{\eta } I(t)= f_{4}(t,M,S,L,I,R)=\tau {L}-(\gamma + \mu +\delta ){I}, \\ {} _{0}^{CF}D_{t}^{\eta }R(t) = f_{5}(t,M,S,L,I,R)=\sigma {L}+\gamma {I}- \mu {R}, \end{cases} $$
(6)

where \(M(0)=\mathbf{N}_{1} \), \(S(0)=S_{0} \), \(L(0)=L_{0} \), \(I(0)=I_{0} \), \(R(0)=R_{0}\). Application of \({}_{0}^{CF}I^{\eta }\) on both sides of (2) gives the following system of integral equations:

$$ \textstyle\begin{cases} M(t)=M(0)+G [f_{1} (t,M,S,L,I,R)-f_{01} ] \\ \hphantom{M(t)=}{} +\overline{G} \int _{0}^{t}f_{1} (\xi ,M,S,L,I,R)\,d\xi , \\ S(t)=S(0)+G [f_{2} (t,M,S,L,I,R)-f_{02} ] \\ \hphantom{S(t)=}{} +\overline{G} \int _{0}^{t}f_{2} (\xi ,M,S,L,I,R)\,d\xi , \\ L(t)=L(0)+G [f_{3} (t,M,S,L,I,R)-f_{03} ] \\ \hphantom{L(t)=}{}+\overline{G} \int _{0}^{t}f_{3} (\xi ,M,S,L,I,R)\,d\xi , \\ I(t)=I(0)+G [f_{4} (t,M,S,L,I,R)-f_{04} ] \\ \hphantom{I(t)=}{}+\overline{G} \int _{0}^{t}f_{4} (\xi ,M,S,L,I,R)\,d\xi , \\ R(t)=R(0)+G [f_{5} (t,M,S,L,I,R)-f_{05} ] \\ \hphantom{R(t)=}{}+\overline{G} \int _{0}^{t}f_{5} (\xi ,M,S,L,I,R)\,d\xi , \end{cases} $$
(7)

where \(G=\frac{(1-\eta )}{ \mathcal{K}(\eta )}\), \(\overline{G}= \frac{\eta }{\mathcal{K}(\eta )} \). Further, we will use \(f_{0i}=f_{i}(0, M(0), S(0), L(0), I(0), R(0))\), \(i=1,2,3,4,5\). Using the initial conditions, we have

$$ \textstyle\begin{cases} M(t)=\mathbf{N}_{1}+G [f_{1} (t,M,S,L,I,R)-f_{01} ] \\ \hphantom{M(t)=}{} +\overline{G} \int _{0}^{t}f_{1} (\xi ,M,S,L,I,R)\,d\xi , \\ S(t)=\mathbf{N}_{2}+G [f_{2} (t,M,S,L,I,R)-f_{02} ] \\ \hphantom{S(t)=}{}+\overline{G} \int _{0}^{t}f_{2} (\xi ,M,S,L,I,R)\,d\xi , \\ L(t)=\mathbf{N}_{3}+G [f_{3} (t,M,S,L,I,R)-f_{03} ] \\ \hphantom{L(t)=}{}+\overline{G} \int _{0}^{t}f_{3} (\xi ,M,S,L,I,R)\,d\xi , \\ I(t)=\mathbf{N}_{4}+G [f_{4} (t,M,S,L,I,R)-f_{04} ] \\ \hphantom{I(t)=}{} +\overline{G} \int _{0}^{t}f_{4} (\xi ,M,S,L,I,R)\,d\xi , \\ R(t)=\mathbf{N}_{5}+G [f_{5} (t,M,S,L,I,R)-f_{05} ] \\ \hphantom{R(t)=}{}+\overline{G} \int _{0}^{t}f_{5} (\xi ,M,S,L,I,R)\,d\xi . \end{cases} $$
(8)

Here, we denote Banach space by \(X=C([0, T]\times \mathscr{R}^{5}, \mathscr{R})\) under the norm

$$ \Vert W \Vert =\max_{t\in [0,T]} \bigl\{ \bigl\vert W(t) \bigr\vert :W=(M,S,L,I,R) \bigr\} , $$

where \(T>0\) such that \(0\leq t\leq T<\infty \).

Theorem 1

(Krasnoselskii fixed point theorem)

Let X be a Banach space and D be a closed and convex subset of X, then there exist two operators A, B for which the following hold:

  1. 1.

    The sum \(Ax+By\) belongs to D;

  2. 2.

    The operator A is a contraction, while the operator B is continuous and compact;

  3. 3.

    at least one solution z in a way that \(Az+Bz=z\) holds.

Let us assume

W(t)= [ M S L I R ] ,W(0)(t)= W 0 = [ N 1 N 2 N 3 N 4 N 5 ] ,F= [ f 1 f 2 f 3 f 4 f 5 ] .

Therefore, system (8) reduces to

$$\begin{aligned}& W(t)=W_{0} +G F(t,M,S,L,I,R)-G F_{0}+\overline{G} \int _{0}^{t} F(\xi ,M,S,L,I,R)\,d\xi, \\& W(t)=W_{0}+G \bigl[F(t,W)-F_{0} \bigr]+ \overline{G} \int _{0}^{t} F(\xi ,W)\,d\xi . \end{aligned}$$
(9)

For further analysis, we suppose that the following assumptions hold:

(\(\mathcal{H}_{1}\)):

There exists \(\mathcal{K}_{F}>0\) for which

$$ \bigl\vert F(t,W)-F(t,\overline{W}) \bigr\vert \leq \mathcal{K}_{F} \vert W-\overline{W} \vert . $$
(\(\mathcal{H}_{2}\)):

There exist \(C_{F}>0\) and \(M_{F}>0\) such that

$$ \bigl\vert F(t,W) \bigr\vert \leq C_{F} \vert W \vert +M_{F}. $$

Theorem 2

In the view of Theorem 1, problem (9) has at least one solution provided \(G \mathcal{K}_{F}\) is less than unity.

Proof

To prove the theorem, we define a compact and closed set D such that \(D={W \in X: \|W\|\leq r}\). Next, we define operators A and B as follows:

$$ \begin{aligned} & AW(t)=W_{0} +G F(t,M,S,L,I,R)-G F_{0}, \\ & BW(t)=\overline{G} \int _{0}^{t} F(\xi ,M,S,L,I,R)\,d\xi . \end{aligned} $$
(10)

To verify that the operator A in (10) is a contraction, we assume \(W,\overline{W} \in X\), so that

$$\begin{aligned} \Vert AW-A\overline{W} \Vert =&\max \bigl\vert AW(t)-A\overline{W} (t) \bigr\vert \\ =&\max \bigl\vert G \bigl[F(t,W)-F(t,\overline{W}) \bigr] \bigr\vert , \end{aligned}$$

from which we have

$$ \Vert AW-A\overline{W} \Vert \leq G K_{F} \bigl[ \Vert W- \overline{W} \Vert \bigr], $$

which clearly indicates that the operator A is a contraction.

Now we show that the operator B is compact and continuous. Consider

$$\begin{aligned} \bigl\vert BW(t) \bigr\vert =& \biggl\vert \overline{G} \int _{0}^{t} F(\xi ,M,S,L,I,R)\,d\xi \biggr\vert \\ \leq &\overline{G} \int _{0}^{t} \bigl\vert F \bigl(\xi ,W(\xi ) \bigr) \bigr\vert \,d\xi . \end{aligned}$$
(11)

Taking max of (11), we have

$$\begin{aligned} \Vert BW \Vert \leq & \overline{G}\max _{t \in [0,T]} \int _{0}^{t} \bigl\vert F \bigl(\xi ,W( \xi ) \bigr) \bigr\vert \,d\xi \\ \leq & \overline{G}\max_{t \in [0,T]} \int _{0}^{t} \bigl[C_{F} \Vert W \Vert +M_{F} \bigr]\,d\xi \\ \leq & \overline{G} T(C_{F} r+M_{F}). \end{aligned}$$
(12)

This implies that B is bounded in (12). Let us assume that in the domain of t we have \(t_{1}< t_{2}\). One may write

$$\begin{aligned} \bigl\vert BW(t_{2} )-BW(t_{1}) \bigr\vert =& \biggl\vert \overline{G} \int _{0}^{t_{2}}F(\xi ,W)\,d\xi -\overline{G} \int _{0}^{t_{1}}F(\xi ,W)\,d\xi \biggr\vert \\ \leq & \overline{G}(C_{F} r+M_{F}) (t_{2}-t_{1}). \end{aligned}$$
(13)

If \(t_{2}\) approaches \(t_{1}\), then the right-hand side of (13) goes to zero. Consequently, \(t_{2} \rightarrow t_{1}\), which leads to

$$ \bigl\vert BW(t_{2} )-BW(t_{1}) \bigr\vert \rightarrow 0. $$

It follows that B is equicontinuous and, consequently, B is compact continuous. This implies that B is a completely continuous operator. Hence, all the conditions of Theorem 1 are satisfied. One may immediately conclude that model (2) has at least one solution. □

Theorem 3

There is a unique solution of the model under consideration (2) if the functions \(f_{1}\), \(i=1,2,3,4\), are continuous and \(\overline{G} K_{F}(1+T)<1\).

Proof

To prove the theorem, we define an operator \(P:X\rightarrow X\) such that

$$ PW(t)=W_{0} +G F(t,W)-GF_{0}+\overline{G} \int _{0}^{t_{1}} F(\xi ,W)\,d\xi . $$

Let \(W, \overline{W} \in X\), we may write

$$\begin{aligned} \bigl\Vert P(W)-P(\overline{W}) \bigr\Vert =& \max \bigl\vert {P(W) (t)-P(\overline{W}) (t)} \bigr\vert \\ \leq & \max \bigl\vert G (F(t, W)-F(t,\overline{W}) \bigr\vert + \overline{G} \max | \int _{0}^{t} \bigl\vert (F(t,W)-F(t,P( \overline{W}) \bigr\vert d \xi \\ \leq & \overline{G} K_{F} \Vert W-\overline{W} \Vert +G K_{F} T \Vert W- \overline{W} \Vert . \end{aligned}$$

It follows that

$$\begin{aligned} \Vert W-\overline{W} \Vert \leq \overline{G} K_{F}(1+T) \Vert W-\overline{W} \Vert . \end{aligned}$$
(14)

Consequently, model (2) under investigation has a unique solution. □

4 Construction of general algorithm for the required solution of the considered model

To derive the series type solution for the considered problem, we take \(\mathcal{K}(\eta )=1\) and apply the Laplace transform on both sides of (2). We construct the following algorithm:

$$ \textstyle\begin{cases} \frac{s \mathscr{L}[M(t)]-M(0)}{s+\eta (1-s)}= \mathscr{L} [\theta \rho -(\alpha +\mu ) M ], \\ \frac{s \mathscr{L}[S(t)]-S(0)}{s+\eta (1-s)}= \mathscr{L} [(1- \theta )\rho +\alpha M-\beta SI- \mu S ], \\ \frac{s \mathscr{L}[L(t)]-L(0)}{s+\eta (1-s)}= \mathscr{L} [\beta SI-( \sigma +\tau +\mu )L ], \\ \frac{s \mathscr{L}[I(t)]-I(0)}{s+\eta (1-s)}= \mathscr{L} [\tau L-( \gamma +\mu +\delta )I ], \\ \frac{s \mathscr{L}[R(t)]-R(0)}{s+\eta (1-s)}= \mathscr{L} [\sigma L+ \gamma I-\mu R ]. \end{cases} $$
(15)

After rearranging terms in (15), we have

$$ \textstyle\begin{cases} \mathscr{L} [M(t) ]= \frac{M(0)}{s}+ \frac{s+\eta (1-s)}{s}\mathscr{L} [\theta \rho -(\alpha +\mu ) M ], \\ \mathscr{L} [S(t) ]= \frac{S(0)}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [(1-\theta )\rho +\alpha M-\beta SI- \mu S ], \\ \mathscr{L} [L(t) ]= \frac{L(0)}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [ \beta SI-(\sigma +\tau +\mu )L ], \\ \mathscr{L} [I(t) ]= \frac{I(0)}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [ \tau L-(\gamma +\mu +\delta )I ], \\ \mathscr{L} [R(t) ]= \frac{R(0)}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [ \sigma L+\gamma I-\mu R ]. \end{cases} $$
(16)

Using the initial conditions of system (2), one has

$$ \textstyle\begin{cases} \mathscr{L} [M(t) ]= \frac{\mathbf{N}_{1}}{s}+ \frac{s+\eta (1-s)}{s}\mathscr{L} [\theta \rho -(\alpha +\mu ) M ], \\ \mathscr{L} [S(t) ]= \frac{\mathbf{N}_{2}}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [(1-\theta )\rho +\alpha M-\beta SI- \mu S ], \\ \mathscr{L} [L(t) ]= \frac{\mathbf{N}_{3}}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [\beta SI-(\sigma +\tau +\mu )L ], \\ \mathscr{L} [I(t) ]= \frac{\mathbf{N}_{4}}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [\tau L-(\gamma +\mu +\delta )I ], \\ \mathscr{L} [R(t) ]= \frac{\mathbf{N}_{5}}{s}+\frac{s+\eta (1-s)}{s} \mathscr{L} [\sigma L+\gamma I-\mu R ]. \end{cases} $$
(17)

Let the solution we compute be in the form of an infinite series as follows:

$$\begin{aligned}& M(t)=\sum_{n=0}^{\infty }M_{n} (t),\qquad S(t)=\sum_{n=0}^{\infty }S_{n} (t), \qquad L(t)=\sum_{n=0}^{\infty }L_{n} (t), \\& I(t)=\sum_{n=0}^{\infty }I_{n} (t),\qquad R(t)=\sum_{n=0}^{\infty }R_{n} (t), \end{aligned}$$

and decompose the nonlinear term SI in terms of Adomian polynomial as follows:

$$ S(t)I(t)=\sum_{n=0}^{\infty }A_{n}(t), $$

where \(A_{n}=\frac{1}{\Gamma {(n+1)}}\frac{d^{n}}{d\lambda ^{n}} [ ( \sum_{k=0}^{n} \lambda ^{k} S_{k} ) (\sum_{k=0}^{n} \lambda ^{k} I_{k} ) ]|_{\lambda =0}\)

$$\begin{aligned}& n=0:\quad A_{0}=S_{0}(t)I_{0}(t), \\& n=1:\quad A_{1}=S_{0}(t) I_{1}(t)+S_{1}(t) I_{0}(t), \\& n=2:\quad A_{2}=S_{0}(t) I_{2}(t)+S_{1}(t) I_{1}(t)+S_{2}(t) I_{0}(t), \\& n=3:\quad A_{3}= S_{0}(t) I_{3}(t)+S_{1}(t) I_{2}(t) +S_{2}(t) I_{1}(t)+S_{3}(t) I_{0}(t), \\& n=4:\quad A_{4}=S_{0}(t) I_{4}(t)+S_{1}(t) I_{3}(t)+S_{2}(t) I_{2}(t)+S_{3}(t) I_{1}(t) +S_{4}(t) I_{0}(t), \\& \cdots \\& n=n:\quad A_{n}=S_{0}(t) I_{n}(t)+S_{1}(t) I_{(n-1)}(t)+\cdots +S_{(n-1)}(t) I_{1}(t)+S_{n}(t) I_{0}(t). \end{aligned}$$

In view of these values, model (8) becomes

$$ \textstyle\begin{cases} \mathscr{L} [\sum_{k=0}^{\infty }M_{k} (t) ]= \frac{\mathbf{N}_{1}}{s} +\frac{s+\eta (1-s)}{s}\mathscr{L} [\theta \rho -(\alpha +\mu ) \sum_{k=0}^{\infty }M_{k} (t) ], \\ \mathscr{L} [\sum_{k=0}^{\infty }S_{k} (t) ] \\ \quad = \frac{\mathbf{N}_{2}}{s} +\frac{s+\eta (1-s)}{s}\mathscr{L} [(1-\theta )\rho +\alpha \sum_{k=0}^{\infty }M_{k} (t) -\beta \sum_{k=0}^{\infty }A_{k}- \mu \sum_{k=0}^{\infty }S_{k} (t) ], \\ \mathscr{L} [\sum_{k=0}^{\infty }L_{k} (t) ]= \frac{\mathbf{N}_{3}}{s} +\frac{s+\eta (1-s)}{s}\mathscr{L} [\beta \sum_{k=0}^{\infty }A_{k}-( \sigma +\tau +\mu )\sum_{k=0}^{\infty }L_{k} (t) ], \\ x \mathscr{L} [\sum_{k=0}^{\infty }I_{k} (t) ]= \frac{\mathbf{N}_{4}}{s} +\frac{s+\eta (1-s)}{s}\mathscr{L} [\tau \sum_{k=0}^{\infty }L_{k} (t)-(\gamma +\mu + \delta )\sum_{k=0}^{\infty }I_{k} (t) ], \\ \mathscr{L} [\sum_{k=0}^{\infty }R_{k} (t) ] \\ \quad = \frac{\mathbf{N}_{5}}{s} +\frac{s+\eta (1-s)}{s}\mathscr{L} [\sigma \sum_{k=0}^{\infty }L_{k} (t)+\gamma \sum_{k=0}^{\infty }I_{k} (t)-\mu \sum_{k=0}^{\infty }R_{k} (t) ]. \end{cases} $$
(18)

Now, comparing terms on both sides of (18), we get the following series of problems.

Case 1. When \(n=0 \), we have

$$ \textstyle\begin{cases} \mathscr{L} [M_{0}(t) ]= \frac{\mathbf{N}_{1}}{s}+ \frac{s+\eta (1-s)}{s}\mathscr{L}[\theta \rho ], \\ \mathscr{L} [S_{0}(t) ]=\frac{\mathbf{N}_{2}}{s}+ \frac{s+\eta (1-s)}{s}\mathscr{L} [(1-\theta )\rho ], \\ \mathscr{L} [L_{0}(t) ]=\frac{\mathbf{N}_{3}}{s}, \\ \mathscr{L} [I_{0}(t) ]=\frac{\mathbf{N}_{4}}{s}, \\ \mathscr{L} [R_{0}(t) ]=\frac{\mathbf{N}_{5}}{s}. \end{cases} $$
(19)

Evaluating the inverse Laplace transform, we get

$$ \textstyle\begin{cases} M_{0}(t)=\mathbf{N}_{1}+ \theta \rho [1+\eta (t-1) ], \\ S_{0}(t)= \mathbf{N}_{2}+(1-\theta )\rho [1+\eta (t-1) ], \\ L_{0}(t)=\mathbf{N}_{3}, \\ I_{0}(t)=\mathbf{N}_{4}, \\ R_{0}(t)=\mathbf{N}_{5}. \end{cases} $$
(20)

Case 2. When \(n=1 \), we have

$$ \textstyle\begin{cases} \mathscr{L} [M_{1}(t) ]= \frac{s+\eta (1-s)}{s} \mathscr{L} [-(\alpha +\mu ) M_{0}(t) ], \\ \mathscr{L} [S_{1}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [(\alpha M_{0}(t)- \beta A_{0}(t) - \mu S_{0}(t) ], \\ \mathscr{L} [L_{1}(t) ]=\frac{s+\eta (1-s)}{s}\mathscr{L} [\beta A_{0}(t)-( \sigma +\tau +\mu )L_{0}(t) ], \\ \mathscr{L} [I_{1}(t) ]=\frac{s+\eta (1-s)}{s}\mathscr{L} [\tau L_{0}(t) -(\gamma +\mu +\delta )I_{0}(t) ], \\ \mathscr{L} [R_{1}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [\sigma L_{0}(t)+ \gamma I_{0}(t)-\mu R_{0}(t) ]. \end{cases} $$
(21)

Evaluating the inverse Laplace transform, we get

$$ \textstyle\begin{cases} M_{1}(t)= [-(\alpha +\mu ) M_{0}(t) ] [1+\eta (t-1) ], \\ S_{1}(t)= (\alpha M_{0}(t)-\beta S_{0}(t) I_{0}(t)-\mu S_{0}(t) ) [1+ \eta (t-1) ], \\ L_{1}(t)= (\beta S_{0}(t) I_{0}(t)-( \sigma +\tau +\mu )L_{0}(t) ) [1+ \eta (t-1) ], \\ I_{1}(t)= (\tau L_{0}(t) -(\gamma +\mu +\delta )I_{0}(t) ) [1+\eta (t-1) ], \\ R_{1}(t)= (\sigma L_{0}(t) +\gamma I_{0}(t) )-\mu R_{0}(t)) [1+\eta (t-1) ]. \end{cases} $$
(22)

Case 3. When \(n=2 \), we have

$$ \textstyle\begin{cases} \mathscr{L} [M_{2}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [-(\alpha +\mu ) M_{1}(t) ], \\ \mathscr{L} [S_{2}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [(\alpha M_{1}(t)- \beta S_{1}(t) - \mu S_{1}(t) ], \\ \mathscr{L} [L_{2}(t) ]=\frac{s+\eta (1-s)}{s}\mathscr{L} [\beta S_{1}(t)-( \sigma +\tau +\mu )L_{1}(t) ], \\ \mathscr{L} [I_{2}(t) ]=\frac{s+\eta (1-s)}{s}\mathscr{L} [\tau L_{1}(t) -(\gamma +\mu +\delta )I_{1}(t) ], \\ \mathscr{L} [R_{2}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [\sigma L_{1}(t)+ \gamma I_{1}(t)-\mu R_{1}(t) ]. \end{cases} $$
(23)

Evaluating the inverse Laplace transform, we get

$$ \textstyle\begin{cases} M_{2}(t)=(\alpha +\mu )^{2} M_{0} [1+2 \eta (t-1)+{ \eta }^{2} (\frac{t^{2}}{2!}-2t+1 ) ], \\ S_{2}(t)= [-\alpha (\alpha +\mu )M_{0}-\beta \{S_{0} (\tau L_{0}-( \gamma +\mu +\delta )I_{0} ) \\ \hphantom{S_{2}(t)=}{}+I_{0}(\alpha M_{0} -\beta S_{0} I_{0}-\mu S_{0}) \} -\mu (\alpha M_{0} - \beta S_{0} I_{0}-\mu S_{0}) ] \\ \hphantom{S_{2}(t)=}{}\times [1+2 \eta (t-1)+{\eta }^{2} (\frac{t^{2}}{2!}-2t+1 ) ], \\ L_{2}(t)= [\beta \{S_{0} (\tau L_{0}-(\gamma +\mu +\delta )I_{0} )+I_{0}( \alpha M_{0} -\beta S_{0} I_{0}-\mu S_{0}) \} \\ \hphantom{L_{2}(t)=}{}- \{(\sigma +\tau +\mu ) ( \{\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \} \} ] \\ \hphantom{L_{2}(t)=}{}\times [1+2 \eta (t-1) +{\eta }^{2} (\frac{t^{2}}{2!}-2t+1 ) ], \\ I_{2}(t)= [\tau \{\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \}- \{(\gamma + \mu +\delta ) (\tau L_{0}-(\gamma +\mu +\delta )I_{0} \} ] \\ \hphantom{I_{2}(t)=}{}\times [1+2 \eta (t-1)+{\eta }^{2} (\frac{t^{2}}{2!}-2t+1 ) ], \\ R_{2}(t)= [\tau \{\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \} \\ \hphantom{R_{2}(t)=}{}- \{(\gamma +\mu +\delta ) (\tau L_{0}-\mu (\sigma L_{0}+\gamma I_{0}- \mu R_{0}) \} ] \\ \hphantom{R_{2}(t)=}{}\times [1+2 \eta (t-1)+{\eta }^{2} ( \frac{t^{2}}{2!}-2t+1 ) ]. \end{cases} $$
(24)

Case 4. When \(n=3 \), we have

$$ \textstyle\begin{cases} \mathscr{L} [M_{3}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [-(\alpha +\mu ) M_{2} ], \\ \mathscr{L} [S_{3}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L} [(\alpha M_{2}- \beta P_{2} - \mu S_{2} ], \\ \mathscr{L} [L_{3}(t) ]=\frac{s+\eta (1-s)}{s}\mathscr{L} [\beta P_{2}-( \sigma +\tau +\mu )L_{2} ], \\ \mathscr{L} [I_{3}(t) ]=\frac{s+\eta (1-s)}{s}\mathscr{L} [\tau L_{2} -( \gamma +\mu +\delta )I_{2} ], \\ \mathscr{L} [R_{3}(t) ]=\frac{s+\eta (1-s)}{s} \mathscr{L}[\sigma L_{2}+ \gamma I_{1}-\mu R_{2} ] \end{cases} $$
(25)

and so on. The other terms may similarly be computed.

Evaluating the inverse Laplace transform, we get

$$ \textstyle\begin{cases} M_{3}(t)=-(\alpha +\mu )^{3} M_{0} [1+3 \eta (t-1)+ \eta ^{2} (\frac{t^{2}}{2!}-2t+1 ) \\ \hphantom{M_{3}(t)=}{} +\eta ^{3} ( \frac{t^{3}}{3!}-3\frac{t^{2}}{2!}+3t-1 ) ] \\ S_{3}(t)= [ \alpha (\alpha +\mu )^{2}M_{0}- \beta [ \{S_{0}\tau \{\beta S_{0} I_{0}- (\sigma +\tau +\mu )L_{0} \} \\ \hphantom{S_{3}(t)=}{} - \{(\gamma + \mu + \delta ) (\tau L_{0}-(\gamma +\mu +\delta )I_{0} ) \} ] \\ \hphantom{S_{3}(t)=}{} - [\alpha (\alpha +\mu )M_{0}-\beta \{S_{0} (\tau L_{0}-(\gamma +\mu + \delta )I_{0} ) \\ \hphantom{S_{3}(t)=}{}+I_{0}(\alpha M_{0} -\beta S_{0} I_{0}- \mu S_{0}) \} -\mu (\alpha M_{0}-\beta S_{0} I_{0}-\mu S_{0}) ] I_{0} \\ \hphantom{S_{3}(t)=}{} -\mu [ - \alpha (\alpha +\mu )M_{0}-\beta \{S_{0} (\tau L_{0}-( \gamma +\mu + \delta )I_{0} ) \\ \hphantom{S_{3}(t)=}{} +I_{0}(\alpha M_{0} -\beta S_{0} I_{0}- \mu S_{0}) \}-\mu (\alpha M_{0} - \beta S_{0} I_{0}-\mu S_{0}) ] \\ \hphantom{S_{3}(t)=}{} \times [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3 \frac{t^{2}}{2!}+3t-1 ) ] \\ \hphantom{S_{3}(t)=}{} + [ (\alpha M_{0}-\beta S_{0} I_{0}- \mu S_{0}) (\tau L_{0} -( \gamma +\mu +\delta )I_{0} ) ] \\ \hphantom{S_{3}(t)=}{} \times [1+3 \eta (t-1)+\eta ^{2} (2{t^{2}}-2t+1 )+\eta ^{3} (\frac{t^{3}}{{3}}-2{t^{2}}+3t-1 ) ], \\ L_{3}(t)=\beta [ \{S_{0}\tau \{\beta S_{0} I_{0}- (\sigma +\tau + \mu )L_{0} \} - \{(\gamma +\mu +\delta ) (\tau L_{0}-(\gamma +\mu + \delta )I_{0} ) \} ] \\ \hphantom{L_{3}(t)=}{} + [-\alpha (\alpha +\mu )M_{0}-\beta \{S_{0} (\tau L_{0}-(\gamma +\mu + \delta )I_{0} ) \\ \hphantom{L_{3}(t)=}{}+I_{0}(\alpha M_{0} -\beta S_{0} I_{0}- \mu S_{0}) \} \\ \hphantom{L_{3}(t)=}{} -\mu (\alpha M_{0} -\beta S_{0} I_{0}-\mu S_{0}) I_{0} ]- [(\sigma + \tau +\mu ) [\beta \{S_{0} (\tau L_{0}-(\gamma +\mu +\delta )I_{0} ) \\ \hphantom{L_{3}(t)=}{} +(\alpha M_{0} -\beta S_{0} I_{0}-\mu S_{0}) I_{0} ]-(\sigma +\tau + \mu ) \{\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \} ] \\ \hphantom{L_{3}(t)=}{} \times [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3 \frac{t^{2}}{2!}+3t-1 ) ] \\ \hphantom{L_{3}(t)=}{} +[ (\alpha M_{0}-\beta S_{0} I_{0}-\mu S_{0}) (\tau L_{0} -(\gamma + \mu +\delta )I_{0} ) \\ \hphantom{L_{3}(t)=}{} \times [1+3 \eta (t-1)+\eta ^{2} (2{t^{2}}-2t+1 )+\eta ^{3} ( \frac{t^{3}}{{3}}-2{t^{2}}+3t-1 ) ], \\ I_{3}(t)=[\tau \{\beta S_{0} (\tau L_{0}-(\gamma +\mu +\delta )I_{0} )+ I_{0}(\alpha M_{0} -\beta S_{0} I_{0}- \mu S_{0}) \}\\ \hphantom{I_{3}(t)=}{}-\{(\sigma +\tau + \mu ) \{ \beta S_{0} I_{0} -(\sigma +\tau +\mu )L_{0} \} \\ \hphantom{I_{3}(t)=}{}- [(\gamma +\mu +\delta ) \{\tau (\beta S_{0}I_{0}-( \sigma +\tau +\mu )L_{0} \}-(\gamma +\mu +\delta ) \\ \hphantom{I_{3}(t)=}{}\times (\tau L_{0} -(\gamma +\mu +\delta )I_{0} ) ] \\ \hphantom{I_{3}(t)=}{}\times [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3 \frac{t^{2}}{2!}+3t-1 ) ], \\ R_{3}(t)= \sigma \{\beta S_{0} (\tau L_{0}-(\gamma +\mu +\delta )I_{0} )+ I_{0}(\alpha M_{0} -\beta S_{0} I_{0}- \mu S_{0}) \} \\ \hphantom{R_{3}(t)=}{} -(\sigma +\tau +\mu ) \{\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \} \\ \hphantom{R_{3}(t)=}{}-[( \gamma +\mu +\delta ) \{\tau (\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \} \\ \hphantom{R_{3}(t)=}{} +\gamma [\tau \{\beta S_{0} I_{0}-(\sigma +\tau + \mu )L_{0} \}- \{( \gamma +\mu +\delta ) (\tau L_{0}-(\gamma +\mu +\delta )I_{0} \} \\ \hphantom{R_{3}(t)=}{} -\mu [\tau \{\beta S_{0} I_{0}-(\sigma +\tau +\mu )L_{0} \} \\ \hphantom{R_{3}(t)=}{}- \{(\gamma + \mu +\delta ) (\tau L_{0}-\mu (\sigma L_{0}+\gamma I_{0}- \mu R_{0}) \} ] \\ \hphantom{R_{3}(t)=}{} \times [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3 \frac{t^{2}}{2!}+3t-1 ) ], \end{cases} $$
(26)

and so on. In this way, next terms of the series solution may be computed. Therefore, we get the required solution as follows:

$$ \textstyle\begin{cases} M(t)=M_{0} (t)+M_{1} (t)+M_{2} (t)+M_{3} (t)+\cdots, \\ S(t)=S_{0} (t)+S_{1} (t)+S_{2} (t)+S_{3} (t)+\cdots, \\ L(t)=L_{0} (t)+L_{1} (t)+L_{2} (t)+L_{3} (t)+\cdots, \\ I(t)=I_{0} (t)+I_{1} (t)+I_{2} (t)+I_{3} (t)+\cdots, \\ R(t)=R_{0} (t)+R_{1} (t)+R_{2} (t)+R_{3} (t)+\cdots . \end{cases} $$
(27)

Theorem 4

Let \(\mathscr{X}\) be the Banach space and \(\mathbf{T} : \mathscr{X}\rightarrow \mathscr{X}\) be a contractive nonlinear operator such that, for all \(W, {\bar{W}} \in \mathscr{X}\), \(\|\mathbf{T}(W)-\mathbf{T}({\bar{W}}) \|_{\mathscr{X}}\leq \kappa \|W-{\bar{w}}\|_{\mathscr{X}}\), \(0<\kappa <1\). Using the Banach contraction principle, T has a unique point W such that \(\mathbf{T}{W}={W}\), where \(W=(x,y,z)\). By applying LADM, the series given in (26) can be written as

$$ W_{n}=\mathbf{T}W_{n-1},\qquad W_{n-1}=\sum _{j=1}^{n-1}W_{j},\quad j=1,2,3, \ldots, $$

and let \(W_{0}=W_{0}\in B_{\varepsilon }(W)\), where \(B_{\varepsilon }(W)={\bar{\mathbf{w}}\in \mathscr{X}:\| \bar{\mathbf{w}}-W\|_{\mathscr{X}}<\varepsilon }\), then one has

  1. (i)

    \(W_{n} \in B_{r}(W)\);

  2. (ii)

    \(\lim_{n\rightarrow \infty }W_{n}=W\).

Proof

The proof of the above theorem can be derived in a way similar to that in [39]. □

5 Numerical results and discussion

In this part of the paper, we present numerical results along with illustration regarding the approximate solution of the model under discussion. We take the approximate values for the parameters as given in Table 1. In light of these values, we get the series solution as follows:

$$\begin{aligned} \textstyle\begin{cases} M(t)=90-4.194[1+\eta (t-1)]+0.192411[1+2 \eta (t-1)+{\eta }^{2}( \frac{t^{2}}{2!}-2t+1)] \\ \hphantom{M(t)=}{} -0.008966 [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3\frac{t^{2}}{2!}+3t-1 ) ], \\ S(t)=400-1823.361[1+\eta (t-1)]+8408.2894[1+2 \eta (t-1)+{\eta }^{2}( \frac{t^{2}}{2!}-2t+1)] \\ \hphantom{S(t)=}{} -39200.2682 [1+3 \eta (t-1)+ \eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3\frac{t^{2}}{2!}+3t-1 ) ] \\ \hphantom{S(t)=}{} -355.6365 [1+3 \eta (t-1)+{\eta }^{2}(2{t^{2}}-6t+3)+{\eta }^{3} (\frac{t^{3}}{{3}}-2{t^{2}}+3t-1 ) ], \\ L(t)=100+1811.44[1+\eta (t-1)]-8488.5579[1+2 \eta (t-1)+{\eta }^{2}( \frac{t^{2}}{2!}-2t+1)] \\ \hphantom{L(t)=}{} +39597.60959 [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3\frac{t^{2}}{2!}+3t-1 ) ] \\ \hphantom{L(t)=}{} +355.6365 [1+3 \eta (t-1)+{\eta }^{2}(2{t^{2}}-6t+3) +{\eta }^{3} (\frac{t^{3}}{{3}}-2{t^{2}}+3t-1 ) ], \\ I(t)=50-2.14545[1+\eta (t-1)]+22.6071 [1+2 \eta (t-1)+{\eta }^{2} (\frac{t^{2}}{2!}-2t+1 ) ] \\ \hphantom{I(t)=}{} -106.7888 [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3\frac{t^{2}}{2!}+3t-1 ) ], \\ R(t)=10+4.0454[1+\eta (t-1)]+61.8304 [1+2 \eta (t-1)+{\eta }^{2} (\frac{t^{2}}{2!}-2t+1 ) ] \\ \hphantom{R(t)=}{} +289.93398 [1+3 \eta (t-1)+\eta ^{2} (\frac{t^{2}}{2!}-2t+1 )+\eta ^{3} (\frac{t^{3}}{3!}-3\frac{t^{2}}{2!}+3t-1 ) ]. \end{cases}\displaystyle \end{aligned}$$
(28)

Now we plot the solution up to five terms as given in (28) in Figures 15, corresponding to different fractional orders.

Table 1 Table of parametric values used for simulations of the problem under consideration

It can be observed from Fig. 1 that the immunized population decreases with different fractional orders at different ratio. In the same way the susceptible population is increasing as shown in Fig. 2. The infected and the latently infected population is also increasing as shown in Figs. 3 and 4, respectively. Because the susceptible population is converted to infected or latently infected. Proper cure is applied, then the recovered population will increase as shown in Fig. 5. The process of increase or decrease is initially fastest at lower fractional order and some time after the process is reversed, and the greater is the fractional order the faster is the increasing or decreasing process of population of the respective compartments. It means that fractional order derivatives can express the behavior more globally. The recovered population gradually increases and converges to equilibrium state as time progresses, as shown in Fig. 6.

Figure 2
figure 2

Graphical representation of approximate solutions up to five terms of immunized population \(M(t)\) at different fractional order of the considered model (2)

Figure 3
figure 3

Graphical representation of approximate solutions up to five terms of susceptible population \(S(t)\) at different fractional order of the considered model (2)

Figure 4
figure 4

Graphical representation of approximate solutions up to five terms of latently infected population \(L(t)\) at different fractional order of the considered model (2)

Figure 5
figure 5

Graphical representation of approximate solutions up to five terms of infected population \(I(t)\) at different fractional order of the considered model (2)

Figure 6
figure 6

Graphical representation of approximate solutions up to five terms of recovered population \(R(t)\) at different fractional order of the considered model (2)

6 Conclusion

We have investigated a biological model of TB under Caputo–Fabrizio fractional derivative. We have also established some sufficient results regarding the existence as well as the uniqueness of solution for the considered problem with the help of fixed point theorems. Further we have used a hybrid type method to compute series type solutions for the proposed model. To the best of our knowledge, the said techniques were very rarely used to handle the analytical solutions of FODEs involving nonsingular derivative of Caputo–Fabrizio type in the past. Further the numerical results have been displayed via graphs indicating that the established technique can be used to handle solution of those FODEs involving CFFD efficiently. Further, the mentioned method can be utilized to investigate more nonlinear problems of FODEs involving CFFD.