**2. Preliminaries**

In this segment, some basic definitions, and notations of fuzzy variables are discussed (see References [25–27]).

**Definition 1. (Fuzzy Number)** *A fuzzy number* ψ *is a convex normalized fuzzy set* ψ *of the real line such that:* 

$$\left\{\mu\_{\widehat{\psi}}(\mathfrak{x}) : \mathfrak{R} \to [0,1], \,\forall \, x \in \mathfrak{R}\right\}$$

*where* μψ *is a membership function and is piecewise continuous.*

**Definition 2. (Triangular Fuzzy Number)** *A triangular fuzzy number* ψ *is a convex normalized fuzzy set* ψ *of the real line such that:*

*(a) There exists exactly one x*0 ∈ *with* μψ(*<sup>x</sup>*0) *(x*0 *is called the mean value of* ψ *), where* μψ *is called the membership function of the fuzzy set.*

*(b)* μψ(*x*) *is piecewise continuous.*

> *The membership function* μψ *of a triangular fuzzy number* ψ = (*<sup>a</sup>*1, *b*1, *<sup>c</sup>*1) *is defined as:*

$$\mu\_{\overrightarrow{\psi}}(\mathbf{x}) = \begin{cases} 0, & \mathbf{x} \le a\_1, \\ \frac{\mathbf{x} - a\_1}{\mathbf{b}\_1 - a\_1}, & a\_1 \le \mathbf{x} \le b\_1, \\ \frac{c\_1 - \mathbf{x}}{c\_1 - b\_1}, & b\_1 \le \mathbf{x} \le c\_1, \\ 0, & \mathbf{x} \ge c\_1. \end{cases}$$

**Definition 3. (Single-Parametric Form of Fuzzy Numbers)** *The triangular fuzzy number*ψ = (*<sup>a</sup>*1, *b*1, *<sup>c</sup>*1) *can be characterized by an ordered pair of functions through the* γ*-cut approach* ψ(γ), ψ(γ) = [(*b*1 − *<sup>a</sup>*1)γ + *a*1, −(*<sup>c</sup>*<sup>1</sup> − *<sup>b</sup>*1)γ + *<sup>c</sup>*1]*, where* γ ∈ [0, 1]*. The* γ*-cut form is well-known as the single-parametric form of fuzzy numbers. It is observed that the lower and upper bounds of the fuzzy numbers satisfy the below statements:*

*(i)* ψ(γ) *is a left-bounded nondecreasing continuous function over* [0, 1].


**Definition 4. (Double-Parametric Form of Fuzzy Number)** *Using the single-parametric form, as discussed in Definition 3, we have* ψ = [ψ(γ), ψ(γ)].

*Now we can write this as crisp with DPF as:*

$$\overline{\psi}(\underline{\gamma}, \beta) = \beta(\overline{\psi}(\underline{\gamma}) - \underline{\psi}(\underline{\gamma})) + \underline{\psi}(\underline{\gamma})$$

*where* γ *and* β ∈ [0, 1].

**Definition 5. (Fuzzy Arithmetic)** *For arbitrary fuzzy numbers x* = *<sup>x</sup>*(γ), *<sup>x</sup>*(γ)*, y* = *y*(γ), *y*(γ) *and scalar m, fuzzy arithmetics are well-defined as below:*

*(i) x* = *y if and only if x*(γ) = *y*(γ) *and x*(γ) = *y*(γ).

$$\begin{array}{ll} \mbox{(i)} & \widetilde{\mathbf{x}} + \widetilde{\mathbf{y}} = \lfloor \underline{\mathbf{x}}(\boldsymbol{\gamma}) \rfloor + \underline{\mathbf{y}}(\boldsymbol{\gamma}) , \widetilde{\mathbf{x}}(\boldsymbol{\gamma}) + \widetilde{\mathbf{y}}(\boldsymbol{\gamma}) \rfloor. \\\mbox{(ii)} & \widetilde{\mathbf{x}} \times \widetilde{\mathbf{y}} = \begin{bmatrix} \min\{ \underline{\mathbf{x}}(\boldsymbol{\gamma}) \times \underline{\mathbf{y}}(\boldsymbol{\gamma}) , \underline{\mathbf{x}}(\boldsymbol{\gamma}) \times \overline{\mathbf{y}}(\boldsymbol{\gamma}) \times \underline{\mathbf{y}}(\boldsymbol{\gamma}) \times \underline{\mathbf{y}}(\boldsymbol{\gamma}) , \ \overline{\mathbf{x}}(\boldsymbol{\gamma}) \times \overline{\mathbf{y}}(\boldsymbol{\gamma}) \}. \\\mbox{max}\{ \underline{\mathbf{x}}(\boldsymbol{\gamma}) \times \underline{\mathbf{y}}(\boldsymbol{\gamma}) , \underline{\mathbf{x}}(\boldsymbol{\gamma}) \times \overline{\mathbf{y}}(\boldsymbol{\gamma}) , \overline{\mathbf{x}}(\boldsymbol{\gamma}) \times \underline{\mathbf{y}}(\boldsymbol{\gamma}) , \ \overline{\mathbf{x}}(\boldsymbol{\gamma}) \times \overline{\mathbf{y}}(\boldsymbol{\gamma}) \}. \\\mbox{(iv)} & k\widetilde{\mathbf{x}} = \begin{cases} [k\overline{\mathbf{x}}(\boldsymbol{\gamma}), k\underline{\mathbf{x}}(\boldsymbol{\gamma})] , \ k < 0 \\\ [k\underline{\mathbf{x}}(\boldsymbol{\gamma}), k\overline{\mathbf{x}}(\boldsymbol{\gamma})] , \ k \ge 0 \end{cases}. \end{array}$$

#### **3. Fractional Reduced Di**ff**erential Transform Method**

Let us take an analytic and *k*-times continuously differentiable function ψ(*<sup>x</sup>*, *t*). Assume that ψ(*<sup>x</sup>*, *t*) is denoted as a product of two functions as ψ(*<sup>x</sup>*, *t*) = *<sup>a</sup>*(*x*)*b*(*t*). From Momani and Odibat [28], this function is written as follows

$$\psi(\mathbf{x},t) = \left(\sum\_{m=0}^{\infty} A(m)\mathbf{x}^m\right) \left(\sum\_{n=0}^{\infty} B(n)t^n\right) = \sum\_{m=0}^{\infty} \sum\_{n=0}^{\infty} F(m,n) \mathbf{x}^m t^n \tag{7}$$

where *<sup>F</sup>*(*<sup>m</sup>*, *n*) = *<sup>A</sup>*(*m*)*B*(*n*) is named as the spectrum of ψ(*<sup>x</sup>*, *t*).

**Lemma 1.** *The fractional reduced di*ff*erential transform (FRDT) of an analytic function* ψ(*<sup>x</sup>*, *t*) *is defined as:*

$$\psi\_k(\mathbf{x}) = \frac{1}{\Gamma(1+\alpha k)} [D\_t^{\alpha k} \psi(\mathbf{x}, t)]\_{t=t\_0} \text{ for } k = 0, 1, 2, \dots \tag{8}$$

*The inverse transform of* ψ*k*(*x*) *is well-defined as*

$$\psi(\mathbf{x},t) = \sum\_{k=0}^{\infty} \psi\_k(\mathbf{x}) (t - t\_0)^{ak} \tag{9}$$

*From Equations (8) and (9), we obtain:*

$$\psi\left(\mathbf{x},t\right) = \sum\_{k=0}^{\infty} \frac{1}{\Gamma(\alpha k + 1)} \left[D\_t^{\alpha k} \psi(\mathbf{x},t)\right]\_{t=t\_0} (t - t\_0)^{\alpha k} \tag{10}$$

*In particular, at t*0 = 0*, we have:*

$$\psi(\mathbf{x},t) = \sum\_{k=0}^{\infty} \frac{1}{\Gamma(\alpha k + 1)} [D\_t^{\alpha k} \psi(\mathbf{x},t)]\_{t=0} t^{\alpha k} \tag{11}$$

.

**Theorem 1.** *Let* ψ(*<sup>x</sup>*, *t*), ξ(*<sup>x</sup>*, *t*)*, and* ζ(*<sup>x</sup>*, *t*) *be three analytical functions such that*ψ(*<sup>x</sup>*, *t*) = *R*−<sup>1</sup> *D* [ψ*k*(*x*)]*,* ξ(*<sup>x</sup>*, *t*) = *R*−<sup>1</sup> *D* [ξ*k*(*x*)]*, and* ζ(*<sup>x</sup>*, *t*) = *R*−<sup>1</sup> *D* [ζ*k*(*x*)]*. Hence from References [29–32]):*


*(iii) If* ψ(*<sup>x</sup>*, *t*) = *xmtn, then* ψ*k*(*x*) = *x<sup>m</sup>* δ(*k* − *n*) *where* δ(*k*) = #1, *k* = 0 0, *k* - 0

$$f(i\omega) \quad \text{If } \psi(\mathbf{x}, t) = \mathbf{x}^m t^n \xi(\mathbf{x}, t) \text{, then } \psi\_k(\mathbf{x}) = \mathbf{x}^m \xi\_{k-n}(\mathbf{x}) \dots$$

$$f(\mathbf{z}) \quad \text{If } \boldsymbol{\psi}(\mathbf{x},t) = \boldsymbol{\xi}(\mathbf{x},t)\boldsymbol{\zeta}(\mathbf{x},t) \text{, then } \boldsymbol{\psi}\_k(\mathbf{x}) = \sum\_{i=0}^j \boldsymbol{\xi}\_i(\mathbf{x})\boldsymbol{\zeta}\_{j-i}(\mathbf{x}) = \sum\_{i=0}^j \boldsymbol{\zeta}\_i(\mathbf{x})\boldsymbol{\xi}\_{j-i}(\mathbf{x}).$$

$$\text{(vi)}\quad \text{If } \psi(\mathbf{x},t) = \xi\left(\mathbf{x},t\right)\zeta\left(\mathbf{x},t\right)\zeta\left(\mathbf{x},t\right), \\ \text{then } \psi\_k(\mathbf{x}) = \sum\_{j=0}^k \sum\_{i=0}^j \xi\_i\left(\mathbf{x}\right)\zeta\_{j-i}\left(\mathbf{x}\right)\zeta\_{k-j}\left(\mathbf{x}\right).$$

**Theorem 2.** *Let* ψ(*<sup>x</sup>*, *t*) *and* ξ(*<sup>x</sup>*, *t*) *are two analytical functions such that* ψ(*<sup>x</sup>*, *t*) = *R*−<sup>1</sup> *D* [ψ*k*(*x*)]*, and* ξ(*<sup>x</sup>*, *t*) = *R*−<sup>1</sup> *D* [ξ*k*(*x*)]*. Hence:*

*(i) If* ψ(*<sup>x</sup>*, *t*) = ∂*m* ∂*x<sup>m</sup>* ξ(*<sup>x</sup>*, *t*)*, then* ψ*k*(*x*) = ∂*m* ∂*x<sup>m</sup>* ξ*k*(*x*).

$$\text{(ii)}\quad\text{If }\psi(\mathbf{x},t)=\frac{\partial^{\text{an}}}{\partial t^{\text{an}}}\xi(\mathbf{x},t), \text{ then }\psi\_k(\mathbf{x})=\frac{\Gamma(1+(k+n)\alpha)}{(1+k\alpha)}\xi\_{k+n}(\mathbf{x})\dots$$

**Corollary 1.** *If* ζ(*<sup>x</sup>*, *t*) = *e*δ*t*+<sup>θ</sup> *x, then* ζ*k*(*x*) = δ*k k*!*e*<sup>θ</sup> *x*.

**Corollary 2.** *If* ξ(*<sup>x</sup>*, *t*) = sin(θ*t* + μ*x*) *and* ζ(*<sup>x</sup>*, *t*) = cos(<sup>θ</sup>*<sup>t</sup>* + μ*<sup>x</sup>*)*, then* ξ*k*(*x*) = θ*k k*! sin- *k*π 2 + μ *x and* ζ*k*(*x*) = θ*k k*! cos- *k*π 2 + μ *x* .

In order to explain the concept of FRDTM, let us consider the following equation in the operator form as:

$$L\psi(\mathbf{x},t) + R\psi(\mathbf{x},t) + N\psi(\mathbf{x},t) = h(\mathbf{x},t) \tag{12}$$

with IC:

$$
\psi(\mathbf{x},0) = \mathbf{g}(\mathbf{x})\tag{13}
$$

where *L* = ∂α ∂*t*α ; *R*, *N* are linear, nonlinear operators; and *h*(*<sup>x</sup>*, *t*) is an inhomogeneous source term. Using Theorem 2 and Equations (8) and (12), this reduces to:

$$\frac{\Gamma(1+ak+a)}{\Gamma(1+ak)}\psi\_{k+1}(\mathbf{x}) = H\_k(\mathbf{x}) - R\psi\_k(\mathbf{x}) - N\psi\_k(\mathbf{x}) \text{ for } k = 0, 1, 2\dots \tag{14}$$

where ψ*k*(*x*) and *Hk*(*x*) are the transformed form of ψ(*<sup>x</sup>*, *t*) and *h*(*<sup>x</sup>*, *t*), respectively.

Applying FRDTM on the IC, we obtain:

$$
\psi\_0(\mathbf{x}) = \mathcal{g}(\mathbf{x}) \tag{15}
$$

Using Equations (14) and (15), ψ*k*(*x*) for *k* = 1, 2, 3, ... can be determined. Then,takingthe

$$\text{The inverse transformation of } \|\psi\_k(\mathbf{x})\|\_{k=0}^{\mathcal{U}} \text{ gives the n-term approximate solution as:}$$

$$\underline{\mathbf{n}}$$

$$\psi\_n(\mathbf{x}, t) = \sum\_{k=0}^n \psi\_k(\mathbf{x}) \, t^{ak} \tag{16}$$

Therefore, the analytical result of Equation (12) is written as ψ(*<sup>x</sup>*, *t*) = lim*n*→∞ψ*n*(*<sup>x</sup>*, *t*).

#### **4. Double-Parametric-Based Solution of an Uncertain FDMM Using FRDTM**

To begin with, by applying the single parametric form, the FDMM is changed to an interval-based FDE. At that moment, by applying the DPF, the interval-based FDE is transformed into an FDMM having two parameters that may control the uncertainty. Finally, FRDTM is then applied to solve the corresponding double parametrized FDMM for obtaining the needed solution in terms of intervals/fuzzy variables.

Equations (5) and (6) can now be modified in single-parametric form as:

 *d*α *dt*αψ(*<sup>t</sup>* ; γ), *d*α *dt*αψ(*<sup>t</sup>* ; γ) = −[(0.02γ + *a*1 − 0.02), (−0.02γ + *a*1 + 0.02)]ψ(*t* ; <sup>γ</sup>),ψ(*<sup>t</sup>* ; γ) +[(0.02γ + *b*1 − 0.02), (−0.02γ + *b*1 + 0.02)]ξ(*t* ; γ), ξ(*t* ; γ) '1 − %(0.01γ + δ − 0.01), (−0.01γ + δ + 0.01)&ξ(*t* ; γ), ξ(*t* ; γ)<sup>2</sup>( + [(0.2γ + *c*1 − 0.2), (−0.2γ + *c*1 + 0.2)] *d*α *dt*α ξ(*t* ; γ), *d*α *dt*α ξ(*t* ; γ) = −[(0.02γ + *a*2 − 0.02), (−0.02γ + *a*2 + 0.02)]ξ(*t* ; γ), ξ(*t* ; γ) + [(0.02γ + *b*2 − 0.02), (−0.02γ + *b*2 + 0.02)]ψ(*t* ; <sup>γ</sup>),ψ(*<sup>t</sup>* ; γ) )1 − %(0.01γ + δ − 0.01), (−0.01γ + δ + 0.01)&ψ(*t* ; <sup>γ</sup>),ψ(*<sup>t</sup>* ; γ)2\* + [(0.2γ + *c*2 − 0.2), (−0.2γ + *c*2 + 0.2)]. (17)

with fuzzy ICs:

$$\left[\underline{\psi}(0;\gamma),\overline{\psi}(0;\gamma)\right] = \left[\underline{\xi}(0;\gamma),\overline{\xi}(0;\gamma)\right] = \left[0.1\gamma - 0.1,\ -0.1\gamma + 0.1\right] \tag{18}$$

where Equations (17) and (18) are in interval form. One can find out the solution of this interval equations, but sometimes it is complicated to handle such types of interval equations. Therefore, one may require the double-parametric form to handle this interval computation. Applying double-parametric form to Equations (17) and (18), we have:

$$\begin{cases} \left\| \theta \left( \frac{d^2}{dt^a} \overline{\psi}(t;\gamma) - \frac{d^2}{dt^a} \underline{\psi}(t;\gamma) \right) + \frac{d^2}{dt^a} \underline{\psi}(t;\gamma) \right\| = \left\{ \beta(0.04 - 0.04\gamma) + 0.02\gamma - 0.02 + a\_1 \right\} \\ \left\| \left( \overline{\psi}(t;\gamma) - \underline{\psi}(t;\gamma) \right) + \underline{\psi}(t;\gamma) \right\} + \left\{ \beta(0.04 - 0.04\gamma) + 0.02\gamma + b\_1 - 0.02 \right\} \\ \left\| \left( \overline{\xi}(t;\gamma) - \underline{\xi}(t;\gamma) \right) + \underline{\xi}(t;\gamma) \right\} \Bigg| \begin{aligned} &1 - \left\{ \beta(0.04 - 0.02\gamma) + 0.01\gamma + \delta - 0.01 \right\} \\ &\left\{ \beta(\overline{\xi}(t;\gamma) - \underline{\xi}(t;\gamma)) + 0.01\gamma - 0.01\gamma \right\} \\ &\left\{ \beta(\overline{\xi}(t;\gamma) - \underline{\xi}(t;\gamma)) + \underline{\xi}(t;\gamma) \right\}^2 \end{cases} \tag{19}$$

$$\begin{cases} \beta \left( \frac{d^2}{dt^2} \overline{\xi}(t;\gamma) - \frac{d^2}{dt^2} \underline{\xi}(t;\gamma) \right) + \frac{d^2}{dt^2} \underline{\xi}(t;\gamma) \right] = \{\beta(0.04 - 0.04\gamma) + 0.02\gamma - 0.02 + a\_2\} \\ \beta \left( \overline{\xi}(t;\gamma) - \underline{\xi}(t;\gamma) \right) + \underline{\xi}(t;\gamma) \Big] + \{\beta(0.04 - 0.04\gamma) + 0.02\gamma + b\_2 - 0.02\} \\ \beta \{\overline{\psi}(t;\gamma) - \underline{\psi}(t;\gamma)\} + \underline{\psi}(t;\gamma) \Big) \Bigg[ \begin{array}{c} 1 - \{\beta(0.02 - 0.02\gamma) + 0.01\gamma + \delta - 0.01\} \\ \left\{\beta[\overline{\psi}(t;\gamma) - \underline{\psi}(t;\gamma)] + 0.01\gamma - 0.01\gamma \end{array} \right\} \\ + \{\beta(0.4 - 0.4\gamma) + 0.2\gamma + c\_2 - 0.2\} . \end{cases} \tag{20}$$

with fuzzy ICs:

$$\begin{aligned} \left\{ \beta \Big( \overline{\psi}(0;\gamma) - \underline{\psi}(0;\gamma) \Big) + \underline{\psi}(0;\gamma) \right\} &= \left\{ \beta \Big( \overline{\xi}(0;\gamma) - \underline{\xi}(0;\gamma) \Big) + \underline{\xi}(0;\gamma) \right\} \\ &= \beta (-0.2\gamma + 0.2) + (0.1\gamma - 0.1) \end{aligned} \tag{21}$$

Let us take:

β*d*α *dt*α ξ(*t* ; γ) − *d*α *dt*α ξ(*t* ; γ)+ *d*α *dt*α ξ(*t* ; γ)= *d*α *dt*α ξ(*t* ; γ, β) β- *d*α *dt*αψ(*<sup>t</sup>* ; γ) − *d*α *dt*αψ(*<sup>t</sup>* ; γ) + *d*α *dt*αψ(*<sup>t</sup>* ; γ) = *d*α *dt*αψ(*<sup>t</sup>* ; γ, β) βψ(*t* ; γ) − ψ(*t* ; γ) + ψ(*t* ; γ) = ψ(*t* ; γ, β) βξ(*t* ; γ) − ξ(*t* ; γ) + ξ(*t* ; γ) = ξ(*t* ; γ, β). %β(0.04 − 0.04γ) + 0.02γ − 0.02 + *<sup>a</sup>*1& = *a*1. %β(0.04 − 0.04γ) + 0.02γ − 0.02 + *<sup>a</sup>*2& = *a*2.

$$\begin{aligned} \{\beta(0.04-0.04\gamma)+0.02\gamma+b\_1-0.02\}&=b\_1. \\ \{\beta(0.04-0.04\gamma)+0.02\gamma+b\_2-0.02\}&=\overline{b\_2}. \\ \{\beta(0.4-0.4\gamma)+0.2\gamma+c\_1-0.2\}&=\overline{c\_1}. \\ \{\beta(0.4-0.4\gamma)+0.2\gamma+c\_2-0.2\}&=\overline{c\_2}. \\ \{\beta(0.02-0.02\gamma)+0.01\gamma+\delta-0.01\}&=\overline{\delta}. \\ \{\beta(\overline{\psi}(0:\gamma')-\underline{\psi}(0:\gamma'))+\underline{\psi}(0:\gamma')\}&=\overline{\psi}(0:\gamma,\beta). \end{aligned}$$

and

$$\left\{\beta(\overline{\xi}(0;\gamma) - \underline{\xi}(0;\gamma)) + \underline{\xi}(0;\gamma)\right\} = \overline{\xi}(0;\gamma,\beta)$$

Substituting all the above equations in Equations (19)–(21), we get:

$$\begin{split} \frac{d^a \overline{\psi}(t; \gamma, \emptyset)}{dt^a} &= -\widetilde{a\_1} \widetilde{\psi}(t; \gamma, \emptyset) + \widetilde{b\_1} \widetilde{\xi}(t; \gamma, \emptyset) \Big(1 - \widetilde{\delta \xi^2}(t; \gamma, \emptyset)\Big) + \widetilde{c\_1}. \\\frac{d^a \widetilde{\xi}(t; \gamma, \emptyset)}{dt^a} &= -\widetilde{a\_2} \widetilde{\xi}(t; \gamma, \emptyset) + \widetilde{b\_2} \widetilde{\psi}(t; \gamma, \emptyset) \Big(1 - \widetilde{\delta \psi^2}(t; \gamma, \emptyset)\Big) + \widetilde{c\_2}. \end{split} \tag{22}$$

with ICs:

$$\overline{\psi}(0;\gamma,\beta) = \overline{\xi}(0;\gamma,\beta) = \beta(-0.2\gamma + 0.2) + (0.1\gamma - 0.1) = \eta \tag{23}$$

Solving Equation (22) with the ICs in Equation (23), we have ψ 1(*<sup>t</sup>* ; γ, β) and ψ 2(*<sup>t</sup>* ; γ, β) in terms of γ and β. To find the lower and upper bounds of the solutions in single parametric form, we have to substitute β = 0 and β = 1, respectively. Mathematically these are written as:

$$
\overline{\psi}(t;\gamma,0) = \underline{\psi}(t;\gamma), \\
\overline{\xi}(t;\gamma,0) = \underline{\xi}(t;\gamma) \text{ and } \overline{\psi}(t;\gamma,1) = \overline{\psi}(t;\gamma), \\
\overline{\xi}(t;\gamma,1) = \overline{\xi}(t;\gamma).
$$

Applying FRDTM to both sides of Equation (22), and using Theorems 1 and 2, we have:

$$\begin{split} \frac{\Gamma(1+ka+a)}{\Gamma(1+ka)} \widetilde{\psi}\_{k+1}(\boldsymbol{\gamma},\boldsymbol{\beta}) &= -\widetilde{a}\_{1} \widetilde{\psi}\_{k}(\boldsymbol{\gamma},\boldsymbol{\beta}) + \widetilde{b}\_{1} \widetilde{\xi}\_{k}(\boldsymbol{\gamma},\boldsymbol{\beta}) - \widetilde{b}\_{1} \widetilde{\delta} \left( \sum\_{i=0}^{k} \sum\_{j=0}^{i} \widetilde{\xi}\_{i-j} \widetilde{\xi}\_{j} \widetilde{\xi}\_{k-i} \right) + \widetilde{c}\_{1} \delta(k) \\ \frac{\Gamma(1+ka+a)}{\Gamma(1+ka)} \widetilde{\xi}\_{k+1}(\boldsymbol{\gamma},\boldsymbol{\beta}) &= -\widetilde{a}\_{2} \widetilde{\xi}\_{k}(\boldsymbol{\gamma},\boldsymbol{\beta}) + \widetilde{b}\_{2} \widetilde{\psi}\_{k}(\boldsymbol{\gamma},\boldsymbol{\beta}) - \widetilde{b}\_{2} \widetilde{\delta} \left( \sum\_{i=0}^{k} \sum\_{j=0}^{i} \widetilde{\psi}\_{i-j} \widetilde{\psi}\_{j} \widetilde{\psi}\_{k-i} \right) + \widetilde{c}\_{2} \delta(k) \end{split} \tag{24}$$

Where:

$$\delta(k) = \begin{cases} 1, & k = 0\\ 0, & k \neq 0 \end{cases}$$

Using FRDTM on the IC, we get:

$$
\overline{\psi}\_0(\gamma, \beta) = \overline{\xi}\_0(\gamma, \beta) = \eta \tag{25}
$$

Using Equation (25) in Equation (24), the following values of ψ *k* and ξ *k* for *k* = 1, 2, ... are obtained:

$$\begin{array}{l} \overline{\psi}\_{1} = -\overline{a}\_{1}\eta + \overline{b}\_{1}\eta - \overline{b}\_{1}\overline{\delta}\eta^{3} + \overline{c}\_{1}. \\\overline{\xi}\_{1} = -\overline{a}\_{2}\eta + \overline{b}\_{2}\eta - \overline{b}\_{2}\overline{\delta}\eta^{3} + \overline{c}\_{2}. \end{array} \tag{26}$$

$$
\begin{split}
\widetilde{\psi}\_{2} &= -\widetilde{a}\_{1} \big( -\widetilde{a}\_{1}\eta + \widetilde{b}\_{1}\eta - \widetilde{b}\_{1}\widetilde{\delta}\eta^{3} + \widetilde{c}\_{1} \big) + \widetilde{b}\_{1} \big( -\widetilde{a}\_{2}\eta + \widetilde{b}\_{2}\eta - \widetilde{b}\_{2}\widetilde{\delta}\eta^{3} + \widetilde{c}\_{2} \big) \\ & \qquad - 3\widetilde{b}\_{1}\widetilde{\delta}\eta^{2} \big( -\widetilde{a}\_{2}\eta + \widetilde{b}\_{2}\eta - \widetilde{b}\_{2}\widetilde{\delta}\eta^{3} + \widetilde{c}\_{2} \big) . \\
\widetilde{\xi}\_{2} &= -\widetilde{a}\_{2} \big( -\widetilde{a}\_{2}\eta + \widetilde{b}\_{2}\eta - \widetilde{b}\_{2}\widetilde{\delta}\eta^{3} + \widetilde{c}\_{2} \big) + \widetilde{b}\_{2} \big( -\widetilde{a}\_{1}\eta + \widetilde{b}\_{1}\eta - \widetilde{b}\_{1}\widetilde{\delta}\eta^{3} + \widetilde{c}\_{1} \big) . \\
& \qquad - 3\widetilde{b}\_{1}\widetilde{\delta}\eta^{2} \big( -\widetilde{a}\_{1}\eta + \widetilde{b}\_{1}\eta - \widetilde{b}\_{1}\widetilde{\delta}\eta^{3} + \widetilde{c}\_{1} \big). \end{split} \tag{27}$$

Continuing the above procedure, all the values of ψ*i*, <sup>ξ</sup>*<sup>i</sup>*<sup>∞</sup>*i*=<sup>3</sup> can be calculated. Therefore, according to FRDTM, the n-term solutions of Equation (22) with Equation (23) are written as:

$$\begin{aligned} \overleftarrow{\psi}(t;\mathsf{y},\boldsymbol{\beta}) &= \sum\_{k=0}^{n} \overleftarrow{\psi}\_{k}(\mathsf{y},\boldsymbol{\beta})t^{\mathrm{ak}}.\\ \overleftarrow{\xi}(t;\mathsf{y},\boldsymbol{\beta}) &= \sum\_{k=0}^{n} \overleftarrow{\xi}\_{k}(\mathsf{y},\boldsymbol{\beta})t^{\mathrm{ak}}.\end{aligned} \tag{28}$$

Substituting β = 0 and β = 1, the lower and upper bounds of the solution can be calculated, which, respectively, are as follows:

$$\begin{aligned} \widetilde{\psi}(t;\mathcal{V},0) &= \sum\_{k=0}^{n} \widetilde{\psi}\_{k}(\mathcal{V},0)t^{ak}. \\ \widetilde{\xi}(t;\mathcal{V},0) &= \sum\_{k=0}^{n} \widetilde{\xi}\_{k}(\mathcal{V},0)t^{ak}. \end{aligned} \tag{29}$$

and

$$\begin{aligned} \widetilde{\psi}(t;\boldsymbol{\gamma},1) &= \sum\_{k=0}^{n} \widetilde{\psi}\_{k}(\boldsymbol{\gamma},1)t^{ak}. \\ \widetilde{\xi}(t;\boldsymbol{\gamma},1) &= \sum\_{k=0}^{n} \widetilde{\xi}\_{k}(\boldsymbol{\gamma},1)t^{ak}. \end{aligned} \tag{30}$$

#### **5. Results and Discussion**

In this section, an approximate solution of a fuzzy FDMM using FRDTM has been studied. Various numerical computations have been carried out by taking different values of parameters involved in the equation and ICs. In this article, all the figures and tables are included by considering the values of the parameters as *a*1 = 0.05, *b*1 = 0.04, *c*1 = 0.2, *a*2 = 0.07, *b*2 = 0.06, *c*2 = 0.3 and δ = 0.01 (see References [16,17]). The achieved outcomes are compared with the solution of Singh et al. [16] and Goyal et al. [17], which show the validation of the present study. Calculated results are displayed in terms of plots.

Here, all the numerical calculations have been computed by truncating the infinite series to a finite number of terms (*n* = <sup>5</sup>). Fuzzy solutions of FDMM are portrayed in Figures 1–6 by changing time *t* from 0 to 1 and for different values of α. Next, interval solutions for different values of α have been illustrated in Figures 7–12 by considering γ − *cut* 0.4, 0.8 and 1, and varying time *t* from 0 to 10. From these Figures 7–12, one may see that the line at γ = 1 is the central line, and all other solutions are present on both sides of the γ = 1 line.

**Figure 1.** Lower and upper bounds fuzzy solutions of ψ(*t*) at α = 1.

**Figure 2.** Lower and upper bounds fuzzy solutions of ξ(*t*) at α = 1.

**Figure 3.** Lower and upper bounds fuzzy solutions of ψ(*t*) at α = 0.5.

**Figure 4.** Lower and upper bounds fuzzy solutions of ξ(*t*) at α = 0.5.

**Figure 5.** Lower and upper bounds fuzzy solutions of ψ(*t*) at α = 0.75.

**Figure 6.** Lower and upper bounds fuzzy solutions of ξ(*t*) at α = 0.75.

**Figure 7.** Lower and upper bounds interval solutions of ψ(*t*) at α = 1.

**Figure 8.** Lower and upper bounds interval solutions of ξ(*t*) at α = 1.

**Figure 9.** Lower and upper bounds interval solutions of ψ(*t*) at α = 0.5.

**Figure 10.** Lower and upper bounds interval solutions of ξ(*t*) at α = 0.5.

**Figure 11.** Lower and upper bounds interval solutions of ψ(*t*) at α = 0.75.

**Figure 12.** Lower and upper bounds interval solutions of ξ(*t*) at α = 0.75.

Also, it can be seen that the central line (crisp result, i.e., at γ − *cut* = 1) of Figures 7–12 gradually decreased with a decrease in α. Alternatively, we may say that a decrease in the values of α decreased the adoration of man or woman for his/her partner. From Tables 1–6, it is clear that the lower and upper bounds at different values of α were the same at γ = 1 and the obtained results matched with the solutions of Singh et al. [16] and Goyal et al. [17].


**Table 1.** Fuzzy and crisp solution of ψ(*t*) at α = 1 and γ = 1.

**Table 2.** Fuzzy and crisp solution of ξ(*t*) at α = 1 and γ = 1.




**Table 4.** Fuzzy and crisp solution of ξ(*t*) at α = 0.5.



**Table 5.** Fuzzy and crisp solution of ψ(*t*) at α = 0.75.

**Table 6.** Fuzzy and crisp solution of ξ(*t*) at α = 0.75.

