**3. Numerical Scheme and Stability Analysis**

*3.1. Discretization in the Time Direction*

For the discretization of the GFD, we follow the discretization process as discussed in Refs. [19,32]. We split the time interval [0, *τ*] into M equal parts having step size Δ*t* = <sup>T</sup> M, M ∈ <sup>Z</sup>+, and choose the node points as *tr* <sup>=</sup> *<sup>r</sup>*(Δ*t*), *<sup>r</sup>* <sup>=</sup> 0, 1, 2, ... ,M, with starting point *t*<sup>0</sup> = 0. Assuming *w*(*t*) > 0, *γ* ∈ (0, 1), and *z*(*t*) is a monotonic increasing function on [0, <sup>T</sup> ], such that *<sup>η</sup>* <sup>=</sup> *<sup>z</sup>*(*s*); then, *<sup>s</sup>* <sup>=</sup> *<sup>z</sup>*−1(*η*). The discretization of the GFD of *<sup>v</sup>*(*t*) at node *tr* is given as

$$\begin{split} \boldsymbol{\omega} \circ \mathcal{D}\_t^\gamma \boldsymbol{v}(t\_\tau) &= \frac{[\boldsymbol{w}(t\_\tau)]^{-1}}{\Gamma(1-\gamma)} \sum\_{l=1}^r \int\_{t\_{l-1}}^{t\_l} \frac{[\boldsymbol{w}(\boldsymbol{s})\boldsymbol{v}(\boldsymbol{s})]^\top}{[\boldsymbol{z}(t\_l) - \boldsymbol{z}(\boldsymbol{s})]^\gamma} d\boldsymbol{s}, \\ &= \frac{[\boldsymbol{w}(t\_\tau)]^{-1}}{\Gamma(1-\gamma)} \sum\_{l=1}^r \int\_{\boldsymbol{z}(t\_{l-1})}^{\boldsymbol{z}(t\_l)} \frac{1}{[\boldsymbol{z}(t\_r) - \boldsymbol{\eta}]^\gamma} \cdot \frac{d[\boldsymbol{w}(\boldsymbol{z}^{-1}(\boldsymbol{\eta}))\boldsymbol{v}(\boldsymbol{z}^{-1}(\boldsymbol{\eta}))]}{d\boldsymbol{z}^{-1}(\boldsymbol{\eta})} d\boldsymbol{z}^{-1}(\boldsymbol{\eta}), \\ &= \frac{[\boldsymbol{w}(t\_\tau)]^{-1}}{\Gamma(1-\gamma)} \sum\_{l=1}^r \frac{\boldsymbol{w}(t\_l)\boldsymbol{v}(t\_l) - \boldsymbol{w}(t\_{l-1})\boldsymbol{v}(t\_{l-1})}{\boldsymbol{z}(t\_l) - \boldsymbol{z}(t\_{l-1})} \int\_{\boldsymbol{z}(t\_{l-1})}^{\boldsymbol{z}(t\_l)} \frac{1}{[\boldsymbol{z}(t\_r) - \boldsymbol{\eta}]^\gamma} d\boldsymbol{\eta} + \mathcal{R}\_{r\prime} \\ &= \frac{[\boldsymbol{w}(t\_r)]^{-1}}{\Gamma(2-\gamma)} \sum\_{l=1}^r q\_l [\boldsymbol{w}(t\_l)\boldsymbol{v}(t\_l) - \boldsymbol{w}(t\_{l-1})\boldsymbol{v}(t\_{l-1})] + \mathcal{R}\_{r\prime} \end{split} \tag{19}$$

where

$$q\_l = \frac{[z(t\_r) - z(t\_{l-1})]^{1-\gamma} - [z(t\_r) - z(t\_l)]^{1-\gamma}}{z(t\_l) - z(t\_{l-1})}, \ l = 1, 2, \dots, r,\tag{20}$$

and R*<sup>r</sup>* is the truncation error given by

$$\begin{split} \mathcal{R}\_{r} = \frac{[w(t\_{r})]^{-1}}{\Gamma(1-\gamma)} \sum\_{l=1}^{r} \int\_{z(t\_{l-1})}^{z(t\_{l})} \frac{1}{[z(t\_{l})-\eta]^{\gamma}} \left[ \frac{d[w(z^{-1}(\eta))v(x, z^{-1}(\eta))]}{d\eta} \\ & - \frac{w(t\_{l})v(t\_{l}) - w(t\_{l-1})v(t\_{l-1})}{z(t\_{l})-z(t\_{l-1})} \right] d\eta. \end{split} \tag{21}$$

**Lemma 1.** *The coefficient ql*, *l* = 1, 2, ... ,*r*, *given by Equation (20), satisfies qr* > *qr*−<sup>1</sup> > ... > *q*<sup>0</sup> > 0.

**Proof.** For the proof of Lemma 1, see Ref [32].
