Journal of Nonlinear Science (2025) 35:100 https://doi.org/10.1007/s00332-025-10199-8 Adaptive Multidomain Numerical Solution for Singularly Perturbed Fractional Differential Equation: Chebyshev Pseudospectral Method Yusuf Olatunji Tijani1 ·Olumuyiwa Otegbeye2 · Shina Daniel Oloniiju1 Received: 24 October 2024 / Accepted: 28 July 2025 © The Author(s) 2025 Abstract Singularly perturbed differential equations pose significant challenges for many numerical and semi-analytical solution methods in classical calculus. The complexity of these equations increases when dealing with fractional differential equations, which exhibit unique properties, including the nonlocal nature of the arbitrary non-integer order differential operators. This study presents an adaptive multidomain Chebyshev pseudospectralmethod to effectively approximate the solutions of singularly perturbed fractional differential equations. In each subdomain, the fractional differential operator in the Caputo sense accounts for both local derivative effects and memory by dividing the discrete fractional operator into two distinct components, one representing local behavior and the other capturing memory effects. The adaptivity of the method is con- trolled by ensuring that the maximum residual error in each domain remains below a predetermined tolerance value. If the maximum error exceeds the set tolerance, the size of the subinterval is reduced by a specified ratio, and the solution is recomputed within that interval. To assess the performance of the adaptive method, a compara- tive study with the multidomain pseudospectral method with uniform step length is used. The accuracy of the adaptive multidomain Chebyshev pseudospectral method is validated through a series of numerical tests on various fractional differential equa- tions, thus demonstrating the effectiveness of the approach. The adaptive multidomain pseudospectral method offers a robust solution for handling the challenges posed by Communicated by Anthony Bloch. B Yusuf Olatunji Tijani olodyusuf@gmail.com Olumuyiwa Otegbeye olumuyiwa.otegbeye@wits.ac.za Shina Daniel Oloniiju s.oloniiju@ru.ac.za 1 Department of Mathematics, Rhodes University, Makhanda, PO Box 94, Grahamstown 6140, South Africa 2 School of Computer Science and Applied Mathematics, University of the Witwatersrand, Private Bag 3, Braamfontein, Johannesburg 2050, South Africa 0123456789().: V,-vol 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00332-025-10199-8&domain=pdf 100 Page 2 of 22 Journal of Nonlinear Science (2025) 35:100 singularly perturbed differential equations common in fields such aswave propagation, astrophysics, and quantum mechanics. Keywords Adaptive · Singularly perturbed · Multidomain · Fractional differential equation · Spectral method · Errors Mathematics Subject Classification 65L11 · 34A08 · 65L05 · 65L50 1 Introduction Fractional differentiation and integration, which originated around the same time as their classical counterparts, form the foundation for fractional differential equa- tions (Podlubny 1999). In recent years, the understanding and practical application of fractional differential equations has gained significant recognition across various disciplines, including elasticity theory, transport theory, experimental fluid dynamics, and mathematical biology, among others (Ninghu 2020; Diethelm and Freed 1999; Rossikhin and Shitikova 1997). The calculus of arbitrary order has expanded the range of mathematical tools available for modeling real-world phenomena. The oper- ators used in fractional differential equations, like the Riemann–Liouville and Caputo derivatives, are nonlocal, meaning that the derivative at a point depends on all values to the left or right rather than only those in a local neighborhood of the point. This non- local property of fractional operators characterizes fractional differential equations as global in nature (Owolabi and Atangana 2019). A significant advantage of fractional differential equations is their ability to model problems that involve nonlocality and memory effects, aspects that classical differential equations cannot address (Podlubny 1999; Huseynov et al. 2021). It can be argued that most classical differential equa- tions do not readily lend themselves to closed-form solutions (Evans 1998). This may be due to nonlinearity, higher orders, or non-homogeneity of the differential equa- tions. Fractional differential equations are typically more intricate and challenging to solve than classical differential equations (Garrappa 2018). As a result, numerical approximations have emerged as efficient techniques for solving differential equations and addressing uncertainties in differential models. The nonlinearity in the equation itself or in its initial or boundary conditions fundamentally complicates a differential equation (Logan 1994). An example of this complexity is the singularly perturbed differential equation (Kadalbajoo and Patidar 2002). Singularly perturbed problems are a key aspect of initial-boundary value problems and have significant applications in geomechanics, enzyme models, and mathematical biology problems (Kadalbajoo and Patidar 2002). These equations arise when a small parameter multiplies the highest order of the differential equation (Kadalbajoo and Patidar 2002; Endrie and Duressa 2024). Consider the equation − ε d2u(t) dt2 + u(t) = 1, t ∈ (0, 1), u(0) = u(1) = 0. (1) 123 Journal of Nonlinear Science (2025) 35:100 Page 3 of 22 100 FromEq. (1), it can be seen that if the parameter ε is set to zero, it yieldsu(t) = 1,which is inconsistent with the boundary conditions. Consequently, the reduced problem has no solution in C2([0, 1]), the space of twice continuously differentiable functions. This indicates that the solution to Eq. (1) is not well posed when ε becomes very small. This implies that the solution to the problem may exhibit large derivatives near the boundary, and the presence of boundary layers typically makes numerical approxima- tions of singularly perturbed problems challenging (Kadalbajoo and Patidar 2002). It also renders certain classical numerical techniques unstable, resulting in inaccurate numerical results. Finite difference-based methods, such as Crank–Nicholson, may produce reasonable approximations when very fine meshes are used for the grid dis- cretization, but this can lead to a computationally demanding scheme (Sharma et al. 2013). Therefore, developing an efficient, fast, and accurate numerical method for approximating solutions to such differential equations is essential. In this study, we will consider singularly perturbed fractional initial-value differential equations of the form: C D θ 0,t u(t) = g[t, u(t), u(t)2, u′(t), u′′(t), ε]. (2) Here, C Dθ 0,t is the Caputo fractional differential operator of order θ , g : R → R n . Following the study of Shen et al. (2011), spectral methods belong to a class of methods known as the method of weighted residuals, which use infinitely dif- ferentiable functions to approximate the solutions of differential equations. These methods use a truncated series of orthogonal functions, such as Chebyshev, Legen- dre, Jacobi, and Hermite polynomials, to approximate the solutions of differential equations (Canuto et al. 2006). Spectral methods are considered global methods as they discretize differential equations globally to achieve high-order approximations of their solutions (Gottlieb and Gottlieb 2009). Spectral methods can frequently achieve higher accuracy, surpassing the accuracy typically obtained through other traditional numerical methods like the finite difference or finite element methods (Trefethen 2000; Zayernouri and Em Karniadakis 2014). Sweilam and Khader (2010) applied the Legendre pseudospectral method to discretize the fractional advection-dispersion equation involving the Caputo fractional differential operator. Oloniiju et al. (2019) used the Chebyshev pseudospectral method to tackle a continuummechanics problem involving fractional constitutive equations for a fluid dynamics problem. Other appli- cations of the spectral-based methods to classical and fractional differential equations can be found in the following studies (Motsa 2012; Otegbeye 2018; Oloniiju et al. 2021; Doha et al. 2011; Li et al. 2012; Ezz-Eldien and Doha 2023). Ezz-Eldien (2018) applied the spectral tau method, using Jacobi polynomials as basis functions, to solve a system of multi-pantograph equations. The performance of the method was vali- dated through numerical simulations. The spectral Galerkin (SG) technique using the shifted Jacobi polynomials was employed to approximate the solution of a delay time- fractional partial differential equation by Alsuyuti et al. (2023). The effectiveness of the proposed method was assessed by comparing the numerical solutions from five examples with known analytic solutions and results from other techniques in the liter- ature. Spectral-based methods typically use the extrema of orthogonal polynomials as collocation points, which are densely packed near the boundaries. This is particularly 123 100 Page 4 of 22 Journal of Nonlinear Science (2025) 35:100 advantageous for problems like Rayleigh–Bénard convection or channel flow, where a higher density of collocation points is needed to accurately resolve viscous sub- layers (Nnakenyi et al. 2015). However, for problems where rigidity is concentrated in the interior points, the direct application of a single-domain pseudospectral method may lead to some numerical trade-offs. Despite advances in spectral methods, their direct application to certain differential equations, such as fractional singularly per- turbed equations within a single domain, can lead to significant accuracy issues (Motsa et al. 2013). This is due to the nonlinearity introduced by the singular perturbation parameter and the nonlocal nature of these fractional differential equations. To mini- mize the loss of accuracy, domain decomposition or multidomain techniques can be applied (Samuel and Motsa 2019). The multidomain technique for initial-value problems involves dividing the time interval into multiple subintervals and solving the differential equation on each subin- terval sequentially. The process begins with the initial condition, which serves as the starting point or left boundary for the first subinterval. The solution obtained at the last discrete point in the first subinterval is then used as the initial condition for the next subinterval. This procedure continues through all subintervals until the final domain is reached (Oloniiju et al. 2024; Xu and Hesthaven 2014). This approach offers several advantages, such as improving the condition of the algebraic system and allowing larger step sizes compared to single-domain methods. It is particularly effective for solving problems that involve interface tracking and immersed boundary methods (Motsa 2012), and it significantly enhances the accuracy of approximations for stiff and chaotic systems. The multidomain technique is also efficient for problems with long-term dynamics (Oloniiju et al. 2024). However, its direct application to sin- gularly perturbed differential equations may give rise to some significant numerical challenges. To overcome these challenges, an adaptive multidomain approach, which dynamically adjusts each subdomain’s length, may improve efficiency and accuracy. An adaptive numerical technique dynamically adjusts the subdomains based on the complexity of the differential equation or the error estimation of the approximate solu- tion (Peter and Martin 2012). These numerical techniques aim to improve efficiency, accuracy, and computational time by adapting their step sizes, mesh resolutions, or other parameters during the computation (Qureshi et al. 2023). An example of an adaptive numerical technique is adaptive quadrature for numerical integration (Gan- der and Gautschi 2000). Instead of using a fixed number of equally spaced intervals, adaptive quadrature methods automatically refine the intervals in regions where the integrand changes rapidly and coarsen them in smoother areas. This approach enables more accurate results with fewer function evaluations. Another example is adaptive time-stepping ordinary differential equation solvers, like the Runge–Kutta–Fehlberg (RKF) methods, which dynamically adjust their time steps to efficiently handle stiff or non-stiff differential models. This allows the tech- nique to focus computational resources where they are most needed. However, the effectiveness of these methods deteriorates when applied to fractional differential equations. This is because the RKFmethod is local, while fractional differential equa- tions are inherently nonlocal. Global methods, such as the Chebyshev pseudospectral method, are better suited for this purpose (Oloniiju et al. 2024). Adaptive techniques 123 Journal of Nonlinear Science (2025) 35:100 Page 5 of 22 100 have become particularly important for handling nonlinear models, chaotic systems, stiff equations, and singular perturbation problems (Hans-Gorg et al. 1994). A thorough literature review reveals that no existing studies have explored the performance of an adaptive multidomain Chebyshev pseudospectral method for sin- gularly perturbed fractional differential equations. To demonstrate the efficiency of the adaptive multidomain technique, we will use established numerical metrics, such as absolute error, residual error, and computational time, to compare the accuracy and efficiency of this proposed technique with a multidomain pseudospectral method with equal subintervals and other numerical methods. The rest of this manuscript is structured as follows: Section 2 provides the prelim- inaries and definitions required to develop the adaptive multidomain scheme in the context of the calculus of arbitrary orders. Section 3 details the development of the adaptive scheme using the Chebyshev–Gauss–Lobatto quadrature and discusses the development of a nonoverlapping multidomain pseudospectral scheme for differential equations of arbitrary non-integer orders. This section also briefly addresses the error estimates associated with using Chebyshev–Gauss–Lobatto points. Section 4 presents the numerical performance of the method through five distinct numerical examples, and finally, Section 5 concludes the study. 2 Preliminaries and Definitions This section outlines vital preliminary concepts that will be referenced later, estab- lishing the foundation for the terminology used throughout the article. Definition 1 The gamma function, which extends the factorial function, is defined as Podlubny (1999) �(n) = (n − 1)! = ∞∫ 0 tn−1e−tdt, (3) where � is the Euler gamma function. Definition 2 Suppose t > 0 and g ∈ C([0, t]) is continuously defined on a closed interval [0, t]. From Cauchy multiple integrals (Sikora 2023) ∫ t 0 ( ∫ σn 0 ( · · · ( ∫ σ2 0 g(σ1)dσ1 ) · · · ) dσn−1 ) dσn = 1 (n − 1)! ∫ t 0 (t − σ)n−1g(σ )dσ, n ∈ N. (4) Definition 3 (Let θ ∈ R+, t ∈ R+, and g ∈ C([0, t]), then the Riemann–Liouville fractional integral of order, θ is Podlubny (1999)) I θg(t) = 1 �(θ) ∫ t 0 (t − σ)θ−1g(σ )dσ (5) 123 100 Page 6 of 22 Journal of Nonlinear Science (2025) 35:100 Definition 4 Let g(t) ∈ Cn([0, t]), the space of n-times continuously differentiable functions, then the left-sided arbitrary non-integer θ -order differential operator in the Caputo sense is defined as (Oloniiju et al. 2024; Caputo 1967) C D θ 0,t g(t) = ⎧⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ 1 �(n − θ) ∫ t 0 gn(σ ) (t − σ)θ+1−n dσ, for n − 1 < θ < n, dng(t) dtn , for θ = n, n ∈ Z+. (6) InCaputo’s definition of the fractional differential operator, the integer order deriva- tive of the function is computed first, followed by applying the Riemann–Liouville integral of arbitrary order. For θ ∈ (n − 1, n), the Caputo fractional derivative C Dθ 0,t g(t) exists if g ∈ Cn([0, t]). As the Caputo fractional derivative is expressed in integral form, it operates as a nonlocal operator with memory properties, where the current state depends on past states. Remark 1 The Caputo derivative defined in Eq. (6) satisfies several fundamental prop- erties, including linearity, consistency with the classical result for the derivative of a constant, and adherence to both the generalized Leibniz rule and the index law for fractional orders (Podlubny 1999). Definition 5 The fractional derivative of tμ in the Caputo sense is defined as (Li et al. 2012) C D θ 0,t t μ = ⎧⎪⎪⎨ ⎪⎪⎩ 0, μ ∈ N0, μ < �θ� �(μ + 1) �(μ + 1 − θ) tμ−θ , μ ∈ N0, μ ≥ �θ�, (7) where �θ� is the ceiling function of θ and �θ� ≥ θ . 3 Numerical Discretization via Pseudospectral Approximation In this section,wederive themultidomain schemes used in this study.Theprocedure for approximating the solution u(t) of a differential equation of arbitrary real order is first demonstrated, followed by the approximation of the arbitrary real-order derivatives of the function and the subsequent implementation of the multidomain schemes. 3.1 Approximation of the Function u(t) Given the function u(t), which belongs to the space Cn([0,T]), it can be approximated using the shifted Chebyshev polynomials of the first kind. The following recursive expression defines the Chebyshev polynomials: �T,0(t) = 1, 123 Journal of Nonlinear Science (2025) 35:100 Page 7 of 22 100 �T,1(t) = 2t T − 1, �T,2(t) = 8t2 T2 − 8t T + 1, �T,m+1(t) = 2�T,1(t) · �T,m(t) − �T,m−1(t), m ≥ 3. (8) The series expression for the shifted Chebyshev polynomials takes the form: �T,m(t) = m m∑ j=0 (−1)m− j (m + j − 1)!22 j (m − j)!(2 j)!T j t j , m > 0. (9) Hence, using the Chebyshev polynomial, the function u(t) can be approximated in terms of a truncated series as u(t) = Ū(t) � M∑ m=0 �̂m�T,m(t). (10) Here, �̂m are coefficients to be determined and satisfy the orthogonality condition, defined by the following inner product �̂m = 〈Ū(t),�T,m(t)〉 = 1 km ∫ T 0 Ū(t)�T,m(t)√ Tt − t2 dt = k−1 m ∞∑ j=0 Ū(t j )�T,m(t j )λ j , (11) where km = ⎧⎪⎪⎨ ⎪⎪⎩ π m = 0 π 2 otherwise, and λ j = ⎧⎪⎪⎨ ⎪⎪⎩ π 2M j = 0, M π M otherwise. (12) Suppose u(t) is a square-integrable and analytic function; therefore, u(t) is approxi- mated using a linear combination of M + 1 shifted Chebyshev polynomials at M + 1 shifted Chebyshev–Gauss–Lobatto nodes, tr , as u(tr ) � M∑ j=0 Ū(t j )λ j M∑ m=0 1 km �T,m(t j )�T,m(tr ) = M∑ j=0 Ū(t j )Fj (tr ) = FT TŪ, (13) 123 100 Page 8 of 22 Journal of Nonlinear Science (2025) 35:100 where FT = K�T Tλ�T, and λ = ⎡ ⎢⎢⎢⎣ λ0 λ1 . . . λM ⎤ ⎥⎥⎥⎦ , K = ⎡ ⎢⎢⎢⎢⎣ 1 k0 1 k1 . . . 1 kM ⎤ ⎥⎥⎥⎥⎦ , �T = ⎡ ⎢⎢⎢⎣ �T,0(t0) �T,0(t1) . . . �T,0(tM ) �T,1(t0) �T,1(t1) . . . �T,1(tM ) ... ... . . . ... �T,M (t0) �T,M (t1) · · · �T,M (tM ) ⎤ ⎥⎥⎥⎦ . (14) The next step is to find the Caputo fractional derivative of Eq. (13). Assume that u(t) ∈ Cq([0,T]), where M + �θ� + 1 ≤ q, denoting the space of continuously differentiable functions for t > 0 and a fractional order θ > 0, so that C D θ 0,t u(t) ≈ M∑ j=0 Ū(t j )λ j M∑ m=0 1 km �T,m(t j )C D θ 0,t�T,m(t). (15) The Caputo derivative of the space of shifted first kind Chebyshev polynomials, C Dθ 0,t�T,m(t), can be derived using the series expansion of the polynomials (Oloniiju et al. 2024) C D θ 0,t�T,m(t) = m m∑ j=0 (−1)m− j (m + j − 1)!22 j (m − j)!(2 j)!T j C D θ 0,t t j = m m∑ j=�θ� (−1)m− j (m + j − 1)!22 j (m − j)!(2 j)!T j �( j + 1) �( j − θ + 1) t j−θ . (16) The expression, t j−θ , is now written as a linear combination of the shifted Cheby- shev polynomials of the first kind at the Chebyshev–Gauss–Lobatto points. Therefore, Eq. (16) can be written as: C D θ 0,t�T,m(tr ) ≈ D(θ)�T, (17) 123 Journal of Nonlinear Science (2025) 35:100 Page 9 of 22 100 where D(θ) = Z η ZT K. The entries of the matrices Z, η, and are defined by the following: Zm, j = ⎧⎪⎪⎪⎨ ⎪⎪⎪⎩ 1 m = 0, j = 0, m(−1)m− j (m + j − 1)!22m (m − j)!(2 j)!T j m = 1(1)M, j = 0(1)m, 0 otherwise, η j, j = ⎧⎨ ⎩ 0 j < �θ�, �( j + 1) �( j − θ + 1) j = �θ�(1)M, v, j = ⎧⎪⎨ ⎪⎩ 0 j < �θ�,√ πh j−θ+v�( j − θ + v + 1 2 ) �( j − θ + v + 1) j =�θ�(1)M, v=0(1)s, s=0(1)M . ⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ (18) Thus, C Dθ 0,t u(tr ) is given as: C D θ 0,t u(tr ) ≈ (0,TDθ )T Ū, (19) where 0,TDθ = λ �T T KD(θ) �T. The matrix (0,TDθ )T is the fractional differentiation matrixwhich evaluates the fractional derivative of a function at the shiftedChebyshev– Gauss–Lobatto nodes in the domain [0,T]. 3.2 Approximation in Equal Nonoverlapping Partitioned Domains Consider a function u(t) defined on the interval t ∈ [0,T], which is subdivided into p nonoverlapping subintervals. The partitioned points are given by 0 = t0 < t1 < · · · < tp = T, where the length of each subinterval is denoted as δk = tk − tk−1 for k = 1, 2, . . . , p. Assume that t ∈ [tk, tk+1] and θ ∈ (�θ�−1, �θ�], then, the fractional derivative C Dθ 0,t can then be expressed as: C D θ 0,t u(t) = 1 �(�θ� − θ) ∫ t 0 u�θ�(σ ) (t − σ)θ+1−�θ� dσ. (20) This integral representation in Eq. (20) can be further decomposed into a sum of terms corresponding to the ‘memory’ and ‘local’ fractional derivatives in each subinterval as C D θ 0,t u(t) = k∑ f=1 [ C D θ t f −1,t f u(t) ] + C D θ tk ,t u(t), (21) 123 100 Page 10 of 22 Journal of Nonlinear Science (2025) 35:100 where C Dθ t f −1,t f u(t) represents the contribution of the fractional derivative from the previous f -th subinterval and C Dθ tk ,t u(t) is the local fractional derivative in the current subinterval [tk, tk+1]. Eq. (21) is rewritten as C D θ 0,t u(t) = 1 �(�θ� − θ)⎛ ⎝ k∑ f=1 [∫ t f t f −1 u�θ�(σ ) (t − σ)θ+1−�θ� dσ ] + ∫ t tk u�θ�(σ ) (t − σ)θ+1−�θ� dσ ⎞ ⎠ . (22) The first term in Eq. (22) represents the ‘memory’ term of the fractional derivative, incorporating contributions from all previous subintervals. In contrast, the second term is the local fractional derivative within the current subinterval, [tk, tk+1]. These terms are then approximated using the pseudospectral method at M + 1 Chebyshev–Gauss– Lobatto (CGL) points tkr in the subinterval [tk, tk+1]. The local fractional derivative in the subinterval [tk, tk+1] is approximated as: C D θ tk ,t u(tkr ) ≈ (tk ,tk+δkD θ )T Ūk, t ∈ [tk, tk+1], (23) where Ūk is u(t) evaluated at the CGL points within the k-th subinterval. The memory term, C Dθ t f −1,t f u(t), is also approximated as C D θ t f −1,t f u(t) ≈ M∑ j=0 Ū(t fj )λ j M∑ m=0 1 km �δk ,m(t fj ) C D θ t f −1,t f �δk ,m(t), t fj ∈ [t f −1, t f ]. (24) In Eq. (24), the Chebyshev polynomials,�δk ,m(t), are evaluated at theM+1 shifted CGL points, tkr ∈ [tk, tk+1]. Therefore, the definite integrals, C Dθ t f −1,t f �δk ,m(tkr ), defined as Qm,r = ⎧⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎩ 0 j < �θ� m! �(�θ� − θ)(m − �θ�)! ∫ t f t f −1 (σ − t f−1) m−�θ� (tkr − σ)θ+1−�θ� dσ, tkr ∈ [tk, tk+1], j = �θ�, . . . , M, (25) can be evaluated using any quadrature formula. Hence, the approximation of the mem- ory term evaluated at the collocation points, t fr ∈ [t f −1, t f ], becomes k∑ f =1 C D θ t f −1,t f u(t fr ) ≈ k∑ f=1 (t f −1,t f D θ )T Ū f , (26) where t f −1,t f D θ = λ�T δk KZQ and the entries of the matrix,Q, are defined by Eq. (25). 123 Journal of Nonlinear Science (2025) 35:100 Page 11 of 22 100 3.3 Adaptive Nonoverlapping Partitioned Domain for Initial-Value Problem Consider the function, u(t), defined on the interval t ∈ [0,T]. In contrast to the multidomain method with equal subinterval lengths described in Sect. 3.2, where the interval [0,T] is divided into p nonoverlapping subintervals, the adaptive approach dynamically adjusts both the length of each subinterval and the number of subintervals p based on a predefined tolerance and how well the numerical scheme approxi- mates the solution of the differential equation. The subdomains and partition points, 0 = t0, t1, . . . , tp = T, are not predetermined but are generated iteratively as the computation progresses. Beginning with t0 = 0, the length of each subinterval δk = tk − tk−1 is adaptively determined according to the following procedure: 1. Initializing the computation: The computation starts at the initial point t0 = 0, with an initial step length δ1 defined. The solution is computed in the first subinterval [t0, t1]. A tolerance value, denoted as tol, is set to check and control the accuracy of the solution within each subinterval. 2. Estimating the residual error: The norm of the residual error in each subinterval [tk−1, tk] is estimated and checked against the preset tolerance using the following formula: Residual error= ∥∥∥C Dθ tk−1,t u(t)−g[t, u(t), u(t)2, u′(t), u′′(t), ε] ∥∥∥∞ < tol. (27) 3. Adjusting the step length: • If the norm of the residual error is less than or equal to the preset tolerance, the approximate solution in the current subinterval [tk−1, tk] is accepted, and the procedure doubles the latest step length and moves forward to approximate the solution in the next subinterval, [tk, tk+1]. • If the norm of the residual error exceeds the set tolerance, the step length, δk , is reduced using the expression δnewk = 0.9 × δoldk × ( tol Residual error ) 1 M , (28) and the solution is recomputed in the adjusted subinterval. 4. Extending the adaptive domain: The procedure adaptively partitions the interval [0,T] by sequentially creating new subintervals [tk, tk+1] based on the adaptive criterion. The number of subintervals, p, varies as needed, and the process continues until the endpoint tp = T of the domain is reached. To ensure that the partitioning does not exceed the domain of the differential equation, we set the condition that if tk + δk > T, then, δk = T − tk . 123 100 Page 12 of 22 Journal of Nonlinear Science (2025) 35:100 3.4 Error Estimates This subsection presents the general error estimates of the approximation in terms of the Chebyshev polynomials. Assume that �T,mu is the projection of u(t) onto the space of M + 1 Chebyshev polynomials. The error estimate for u − �T,mu in the Sobolev space equipped with its usual norm for u(t) ∈ Hn(0,T) and n ∈ [1,∞) is defined as ‖u − �T,mu‖ 2(0,T) ≤ CM−n|u|Hn,M (0,T). (29) We have the following error estimate for interpolation on the Gauss–Lobatto nodes ‖u − �T,mu‖ 2(0,T) ≤ CM1−n|u|Hn,M (0,T), (30) where C is independent of M and n. The proof of the above estimates is provided in Canuto et al. (2006). 4 Numerical Experimentation and Results This section presents a detailed analysis of the performance of the proposed adaptive multidomain method. Various numerical examples are used to illustrate the accuracy, efficiency, and stability of the technique when applied to singularly perturbed frac- tional differential equations. The results are compared to the nonadaptive multidomain scheme using key metrics such as the absolute error (the difference between the exact and approximate solutions), the norm of the residual error, and the computational time to demonstrate the effectiveness of the adaptive technique. Example 4.1 Consider the singularly perturbed linear fractional differential equation C 0 D θ t u(t) + εu′′(t) + u′(t) = 2 �(3 − θ) t2−θ + 2ε + 2t, for 0 < θ < 1, (31) subject to the initial conditions u(0) = 0, and u′(0) = 0. The closed-form solution of Eq. (31) is u(t) = t2. The right-hand side belongs to the space of at least second-order continuously differentiable functions, i.e., C2([0, T ]). Table 1 compares the maximum absolute errors, residual errors, and computational times between the adaptive multidomain pseudospectral method and the multidomain pseudospectral method with uniform step sizes (i.e., p = 50, 200 and 300) (Oloniiju et al. 2024). The results indicate that the adaptive Chebyshev pseudospectral method yields lower absolute and residual errors than themultidomain technique with uniform step size. This improvement demonstrates the effectiveness of the adaptive approach in achieving higher accuracy. For this comparative analysis, the tolerance of the residual error for the adaptive technique was set at 10−3. The table also suggests that the absolute errors of the equally spaced multidomain method decrease as the number of subdomains increases. However, this improvement comes at the cost of a higher 123 Journal of Nonlinear Science (2025) 35:100 Page 13 of 22 100 Table 1 The ∞-norms of the absolute and residual errors and CPU time of Example 4.1 with M = 10 Metrics ε = 0.01, θ = 0.5, T = 5.0, ξ = 5.0 Adaptive scheme (p = 1) p = 50 p = 200 p = 300 Absolute error 9.3478 × 10−4 2.2464 × 10−2 2.1780 × 10−2 1.9003 × 10−2 Residual error 1.3855 × 10−13 5.6626 × 10−11 2.4131 × 10−10 7.3949 × 10−10 CPU time 0.012s 7.818s 119.515s 269.316s Table 2 The ∞-norms of the absolute and residual errors of Example 4.1 for different values of the perturbation parameter, ε with θ = 0.5, M = 10, and T = 5 Metrics θ = 0.5, ξ = 5.0, T = 5.0, p = 1 ε = 0.1 ε = 0.01 ε = 0.001 ε = 0.0001 Absolute error 1.4324 × 10−3 9.3478 × 10−4 8.9589 × 10−4 8.9209 × 10−4 Residual error 3.5171 × 10−13 1.3855 × 10−13 6.7501 × 10−14 5.0626 × 10−14 Table 3 The maximum absolute and residual errors of Example 4.1 for different initial step size, ξ with ε = 0.01, M = 10, θ = 0.5, and T = 5 Metrics θ = 0.5, ε = 0.01, T = 5.0 ξ = 1.0 ξ = 2.0 ξ = 3.0 ξ = 4.0 Absolute error 8.9307 × 10−3 8.3432 × 10−3 8.1970 × 10−3 5.5004 × 10−3 Residual error 2.2515 × 10−13 1.4521 × 10−13 3.1796 × 10−13 1.0897 × 10−12 Number of subdomains (p) 3 2 2 2 computational time. Tables 2 and 3 present the maximum errors obtained by varying the singularly perturbed parameter ε and the step size of the first domain, ξ . The adaptive method uses only a single domain for the results presented in Table 2. The results in Tables 2 and 3 , respectively, clearly show that the absolute errors improve as ε → 0 and the initial step size increases, i.e., ξ → T . The residual error data in Table 2 highlight that the method proposed approximates the solution of Eq. (31) more accurately as ε approaches zero. Table 3 shows that the adaptive method uses three subdomains when the initial step size equals one. However, when the initial step sizes equal 2, 3, and 4, the method only requires two subdomains. This demonstrates how the adaptive approach adjusts the number of subdomains based on the complexity of the solution in each interval. Example 4.2 Consider the nonlinear singularly perturbed fractional differential equa- tion C 0 D θ t u(t) + εu′′(t) + u′(t) + u2(t) = 6 �(4 − θ) t3−θ − 2 �(3 − θ) t2−θ + ε(6t − 2) + 3t2 − 2t + (t3 − t2)2. (32) 123 100 Page 14 of 22 Journal of Nonlinear Science (2025) 35:100 Table 4 The absolute and residual error norms and CPU time of Example 4.2 with M = 10 Metrics ε = 0.01, θ = 0.5, T = 5.0, ξ = 5.0 Adaptive domain (p = 1) p = 50 p = 200 p = 300 Absolute error 8.5023 × 10−4 1.7049 × 10−2 2.0310 × 10−3 1.7557 × 10−3 Residual error 6.1299 × 10−10 1.2015 × 10−8 1.2119 × 10−9 1.1895 × 10−9 CPU time 0.035s 7.863s 229.459s 536.195s Table 5 The norms of the absolute and residual errors of Example 4.2 using different values of the pertur- bation parameter, ε, with θ = 0.5, M = 10 and T = 5 Metrics θ = 0.5, ξ = 5.0, T = 5.0, p = 1 ε = 0.1 ε = 0.01 ε = 0.001 ε = 0.0001 Absolute error 1.1346 × 10−3 8.5023 × 10−4 8.2330 × 10−4 8.2065 × 10−4 Residual error 6.2937 × 10−10 6.1299 × 10−10 6.0936 × 10−10 6.0254 × 10−10 Table 6 The maximum absolute and residual errors of Example 4.2 for different initial domain length, ξ , with ε = 0.01, M = 10, and T = 5 Metrics θ = 0.5, ε = 0.01, T = 5.0 ξ = 1.0 ξ = 2.0 ξ = 3.0 ξ = 4.0 Absolute error 7.7770 × 10−4 5.0100 × 10−4 4.4339 × 10−4 6.1360 × 10−4 Residual error 1.9707 × 10−8 1.2778 × 10−8 1.9955 × 10−8 1.8087 × 10−8 Number of subdomains (p) 3 2 2 2 Here, θ ∈ (0, 1] with the starting conditions u(0) = 0, and u′(0) = 0. The analytical solution to Eq. (32) subject to its initial conditions is given as u(t) = t2(t − 1). The right-hand side belongs to the space of second-order continuously differentiable functions. The residual and absolute error norms presented in Table 4 mirror the trends observed for Example 4.1. The residual errors are nearly zero, indicating that the numerical scheme closely approximates the solution of the differential equation, albeit with varying accuracy. As shown in Table 4, the adaptive method surpasses the mul- tidomain scheme with uniform step size in terms of error reduction and computational efficiency. The tolerance set for the adaptive scheme for this example is 10−3. A similar trend to that observed in Table 2 is evident in Table 5, where the pro- posed method effectively solves the arbitrary real-order differential equation as the perturbation parameter approaches zero. However, unlike the consistent reduction in absolute error with increase in domain length seen in Table 3, Table 6 reveals a distinct pattern. Specifically, the absolute errors decrease for ξ = 1.0 to ξ = 3.0, but subse- quently increase as ξ progresses from 3.0 to 5.0. When the initial domain length is 1.0 (i.e., ξ = 1.0), the adaptive method uses three subdomains (p = 3). By contrast, for ξ = 2.0, ξ = 3.0, and ξ = 4.0, only two subdomains (p = 2) are required. 123 Journal of Nonlinear Science (2025) 35:100 Page 15 of 22 100 Table 7 The maximum absolute and residual errors and CPU time for Example 4.3 with M = 10 Metrics ε = 0.01, θ = 0.5, T = 5.0, ξ = 5.0 Adaptive domain (p = 1) p = 50 p = 200 p = 300 Absolute error 3.4525 × 10−5 4.7580 × 101 1.2889 × 103 1.1311 × 102 Residual error 2.4364 × 10−11 1.9016 × 102 2.2390 × 10−3 1.6084 × 10−8 CPU time 0.014s 8.405s 124.019s 281.111s Table 8 The maximum absolute and residual errors of Example 4.3 for different values of the perturbation parameter, ε, with θ = 0.5, M = 10 and T = 5 Metrics θ = 0.5, ξ = 5.0, T = 5.0, p = 1 ε = 0.1 ε = 0.01 ε = 0.001 ε = 0.0001 Absolute error 3.4525 × 10−4 3.4525 × 10−5 3.4525 × 10−6 3.4525 × 10−7 Residual error 2.4364 × 10−9 2.4364 × 10−11 2.4364 × 10−13 2.4364 × 10−15 Example 4.3 Given the nonlinear singularly perturbed fractional differential equation ε C 0 D θ t u(t) + u2(t) = ε2 �(2 − θ) t1−θ + ε2(1 + t)2, for 0 < θ ≤ 1, (33) where u(0) = ε. The exact solution to Eq. (33) is given as u(t) = ε(1 + t). A detailed analysis of Table 7 underscores the importance of using the adaptive multidomain technique. Notably, only the solution from the adaptive multidomain Chebyshev pseudospectral method agrees with the closed-form solution with an accu- racy of 10−5. In contrast, the multidomain pseudospectral method with uniform step size is unable to accurately resolve the singularly perturbed fractional differential equation in Example 4.3, resulting in an absolute error of 1.1311× 102 with 300 sub- domains. Additionally, Table 8 confirms that the performance of the adaptive method improves as the value of the perturbation parameter approaches zero. The tolerance set for the results in Table 9 is 10−1; the method fails to converge for smaller toler- ances. Furthermore, for this particular example, using a multidomain approach does not yield an accurate solution to the differential equation, as evidenced by the results in Tables 7 and 9 . These tables illustrate that while effective in other cases, the multidomain technique struggles to provide sufficient accuracy for this problem. The inaccuracies, even when using multiple subdomains, indicate that the structure of this differential equation requires a numerical scheme in a single domain (Table 8), as a multidomain technique, in this instance, does not improve the solution quality. Example 4.4 Consider the singularly perturbed fractional differential equation, which is an interior layer problem (Ramos 2005) ε C 0 D θ t u(t) + u(t) = g(t) + εg′(t), for 0 < t, θ ≤ 1, (34) 123 100 Page 16 of 22 Journal of Nonlinear Science (2025) 35:100 Table 9 The norms of the absolute and residual errors of Example 4.3 for different starting domain length, ξ , with ε = 0.01 and T = 5 Metrics θ = 0.5, ε = 0.01, T = 5.0 ξ = 1.0 ξ = 2.0 ξ = 3.0 ξ = 4.0 Absolute error 3.9980 × 10−1 1.7007 × 10−1 1.6734 × 10−1 7.0310 × 10−1 Residual error 6.6713 × 10−2 5.3586 × 10−2 3.2314 × 10−2 6.9536 × 10−1 Number of subdomains (p) 3 2 2 2 where g(t) = 10 − (10 + t)e−t and u(0) = 10. Typically, this type of differential equation is characterized by solutions with steep gradients or layers that become increasingly thin as the perturbation parameter decreases. When the perturbation parameter ε is small, the solution u(t) exhibits a boundary layer near t = 0 or another critical point, where the solution changes rapidly. Outside of this boundary layer, the solution changes more gradually. The challenge in numerical methods for this kind of problem is accurately capturing the dynamics in the boundary layer while maintaining the overall accuracy of the solution. The analytical solution of Eq. (34), for the case when θ = 1.0, is given as u(t) = 10 − (10 + t)e− t ε + 10e− t ε . Figure 1a illustrates the dynamics of the solution of Eq. (34) as ε → 0, highlighting the presence of a steep gradient t = 0. As ε decreases beyond 10−3, the solutions of the equation begin to superimpose, and no significant differences are observed. Figure 1b depicts the norms of the absolute and residual errors in each subdomain for different values of the perturbation parameter. The number of subdomains used remains consistent when the value of the perturbation parameter is changed, while the initial step length remains constant. For instance, when ξ = 0.01, seven subdomains were used for the different values of ε. According to Table 10, for larger ε, the choice of initial domain length ξ has less impact on the accuracy of the approximate solution. In contrast, smaller ε generally benefits from a smaller initial domain length ξ to achieve lower errors. The table also shows that the accuracy of the approximate solution deteriorates as ε → 0. Table 11 presents the maximum residual error for two different non-integer orders of the differential equation, noting that the residual errors were computed as the exact solutions for these cases are unknown. Table 11 shows that for both values of θ , decreasing ξ generally leads tomore signif- icant residual errors,while smaller values of ε result inmoreminor residual errors. This suggests that themethod exhibits greater accuracy for smaller perturbation parameters. As θ → 1, the residual errors decrease, which indicates that solving the fractional differential equation is more challenging than the classical case. Additionally, it was observed that for a given initial domain length, the number of subdomains used by the adaptive method remains constant, regardless of changes in the perturbation parameter for Example 4.4. The adaptivity criterion was only met once across all the subdomains used in this example. 123 Journal of Nonlinear Science (2025) 35:100 Page 17 of 22 100 Fig. 1 a The exact and numerical solutions of Eq. (34), b the norms of the residual and absolute errors for different values of ε with θ = 1.0, ξ = 10−2, T = 1.0, M = 15 Table 10 The maximum absolute errors for Example 4.4 using different values of the initial domain length ξ and M = 15, θ = 1.0, tol = 10−3 ε M = 15, θ = 1.0, T = 1.0 ξ = 10−1 ξ = 10−2 ξ = 10−3 ξ = 10−4 ξ = 10−5 100 1.698 × 10−8 2.081 × 10−8 2.664 × 10−8 1.880 × 10−8 − 10−1 3.086 × 10−7 9.961 × 10−8 1.036 × 10−7 1.147 × 10−7 8.458 × 10−8 10−2 8.944 × 10−4 3.225 × 10−7 1.437 × 10−7 1.389 × 10−7 1.207 × 10−7 10−3 5.990 × 101 8.746 × 10−4 3.058 × 10−7 1.423 × 10−7 1.503 × 10−7 Number of subdomains (p) 4 7 10 14 17 Table 11 The maximum residual errors for Example 4.4 using different values of the initial domain length ξ , M = 15 and tol = 10−3 ε M = 15, T = 1.0 ξ = 10−1 ξ = 10−2 ξ = 10−3 ξ = 10−4 θ = 0.90 510−1 6.575 × 10−8 8.360 × 10−7 4.025 × 10−6 2.574 × 10−5 10−2 2.424 × 10−9 9.449 × 10−8 1.434 × 10−5 − 10−3 2.072 × 10−10 1.296 × 10−8 2.539 × 10−5 − θ = 0.95 100 3.680 × 10−7 2.609 × 10−6 2.480 × 10−5 1.569 × 10−4 10−1 2.811 × 10−8 3.695 × 10−7 3.541 × 10−6 2.065 × 10−5 10−2 1.780 × 10−9 2.342 × 10−8 5.668 × 10−7 3.567 × 10−5 10−3 7.016 × 10−11 4.308 × 10−9 1.880 × 10−6 1.246 × 10−6 Number of subdomains (p) 4 7 10 14 123 100 Page 18 of 22 Journal of Nonlinear Science (2025) 35:100 Fig. 2 a The plot of the exact and approximate solutions with ξ = 10−1, T = 35.0, and b the number of adaptivity in each subdomain with ε = 10−3, T = 1.0 for θ = 1.0, M = 12 Example 4.5 Next, we evaluate the performance of the adaptive method for the non- linear singularly perturbed fractional Ricatti equation, given as (Ramos 2005) ε C 0 D θ t u(t) + 1 80 u(t) ( u(t) − 20 ) = 0, for 0 < θ ≤ 1, t ∈ R +, (35) with u(0) = 1. The closed-form solution when θ = 1 is u(t) = 20 1 + 19e− t 4ε . Figure 2 demonstrates the accuracy of the adaptive method in capturing the dynam- ics of Example 4.5, where the solution plateaus at u(t) = 20. As the singularly perturbed parameter decreases, the solution reveals an interior layer near the left boundary. In this example, the number of subdomains the adaptive method uses fluc- tuates with changes in ε. Table 12 presents the maximum absolute errors obtained by the adaptive multidomain method when approximating Equation (35) with θ = 1.0. Notably, the error worsens as the perturbation parameter approaches zero. However, reducing the initial step size ξ improves the error, indicating increased accuracy in the approximate solution. In a multidomain method with a uniform step size of ξ = 10−4, the domain [0, 1] would be partitioned into 100, 000 subdomains. In contrast, the adaptive technique resolves the problem using only 42 subdomains, demonstrating its computational efficiency. 5 Conclusion This study examined the performance of an adaptive multidomain pseudospectral method for solving singularly perturbed fractional differential equations. These equa- tions involve derivatives of arbitrary non-integer order with singular perturbations, where small parameters influence the behavior of the solution, often leading to solu- tions with steep gradients. The adaptive method is constructed by projecting the solutionof the differential equationonto the spaceofChebyshevpolynomials andusing 123 Journal of Nonlinear Science (2025) 35:100 Page 19 of 22 100 Table 12 The maximum absolute errors for Example 4.5 using different values of the initial domain length ξ and M = 12, θ = 1.0, tol = 10−3 ε M = 12, θ = 1.0, T = 1.0 ξ = 10−1 ξ = 10−2 ξ = 10−3 ξ = 10−4 ξ = 10−5 100 1.98× 10−11, p = 4 2.19× 10−11, p = 7 3.37× 10−11, p = 10 6.53× 10−11, p = 14 5.37× 10−11, p = 17 10−1 1.73 × 10−9, p = 4 1.84 × 10−9, p = 7 1.05 × 10−9, p = 10 1.15 × 10−9, p = 14 2.64 × 10−9, p = 17 10−2 1.67 × 10−4, p = 4 1.21 × 10−5, p = 7 1.03 × 10−6, p = 10 9.49 × 10−6, p = 14 5.53 × 10−6, p = 17 10−3 2.79 × 10−1, p = 30 1.85 × 10−4, p = 31 1.27 × 10−5, p = 30 9.56 × 10−7, p = 36 1.13 × 10−7, p = 42 the Gauss–Lobatto points to collocate within each domain. The method’s accuracy is assessed through various linear and nonlinear differential equations. The findings of this study reveal the following key points: • The adaptive multidomain scheme demonstrates superior accuracy and computa- tional efficiency performance compared to multidomain with fixed step sizes. • In the case of differential equations with interior layers, the adaptive multidomain scheme performed reasonably well, effectively capturing the solution’s complex behavior while maintaining accuracy across the layers. • As the singularly perturbed parameter approaches zero, the adaptive scheme suc- cessfully handles the presence of interior layers, showcasing its ability to adjust to regions of rapid variation in the solution, ensuring accurate approximation in these regions. • Although the adaptive method generally enhances accuracy and improves com- putational efficiency, the relationship between the initial domain length and the accuracy of the solution remains uncertain, as some problems considered in this study required a large initial step length and others required a smaller step length to achieve the desired accuracy. Further investigation is needed to understand this balance fully. • For the problems considered in this study, using a small tolerance value often results in residual errors that are approximately as small as those obtained with a more significant tolerance. The adaptive multidomain scheme is generally effective for solving singularly perturbed fractional differential equations. In the future, we aim to broaden the applica- bility and improve theperformanceof the proposedmethod through twokeyobjectives: • Refinement of the adaptive strategy:Wewill focus on improving the core algorithm by developing more sophisticated error estimators and mesh-refinement criteria. The goal is to achieve higher levels of accuracywhilemaintaining or even reducing computational cost. • Extension to new problem classes: We plan to adapt the technique to address a wider class of challenging fractional differential equations, notably singularly 123 100 Page 20 of 22 Journal of Nonlinear Science (2025) 35:100 perturbed fractional boundary value problems and fractional partial differential equations. Successfully adapting the method to boundary value problems will require a fun- damental shift from the way the continuity condition is currently enforced across the subdomains. While initial-value problems ensure continuity sequentially from one subdomain to the next, boundary value problems involve handling the subdo- mains globally. This will involve formulating and solving the problem globally to ensure that the boundary conditions at both ends of the domain are satisfied in a single, unified computational step. Author Contributions YOT was involved in conceptualization, initial draft and writing, methodology, soft- ware. OO helped in conceptualization, methodology, review and editing, supervision. SDO contributed to conceptualization, review and editing, supervision, methodology, software. Funding Open access funding provided by Rhodes University. This work is based on the research supported in part by theNational Research Foundation of SouthAfrica (RefNumber TTK2204163593).YusufOlatunji Tijani acknowledges the support of Rhodes University. Data Availability No datasets were generated or analyzed during the current study. Code Availability The code used to generate the results of this study can be made available upon request. Declarations Conflict of interest The authors declare that there are no conflict of interest. The findings and conclusions expressed in this article are solely those of the authors. OpenAccess This article is licensedunder aCreativeCommonsAttribution 4.0 InternationalLicense,which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References Alsuyuti, M.M., Doha, E.H., Bayoumi, B.I., Ezz-Eldien, S.S.: Robust spectral treatment for time-fractional delay partial differential equations. Comput. Appl. Math. 42, 63–73 (2023) Canuto, C., Hussaini,M.Y.,Quarteroni, A., Zang, T.A.: SpectralMethods: Fundamentals in SingleDomains. Springer, Berlin (2006) Caputo, M.: Linear models of dissipation whose q is almost frequency independent-II. Geophys. J. R. Astr. Soc 13, 529–539 (1967) Diethelm, K., Freed, A.D.: On the solution of nonlinear fractional-order differential equations used in the modeling of viscoplasticity. In: Scientific Computing in Chemical Engineering II, pp. 217–224 (1999) Doha, E.H., Bhrawy, A.H., Ezz-Eldien, S.S.: Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations. Appl. Math. Model. 35(12), 5662–5672 (2011) Endrie, N.A., Duressa, G.F.: Numerical method for second order singularly perturbed delay differential equations with fractional order in time via fitted computational method. Partial Differ. Equ. Appl. Math. 10, 100717 (2024) Evans, L.C.: Partial Differential Equations. American Mathematical Society, Providence (1998) 123 http://creativecommons.org/licenses/by/4.0/ Journal of Nonlinear Science (2025) 35:100 Page 21 of 22 100 Ezz-Eldien, S.S.: On solving systems of multi-pantograph equations via spectral tau method. Appl. Math. Comput. 321, 63–73 (2018) Ezz-Eldien, S.S., Doha, E.H.: Fast and precise sepctral method for solving pantograph type Volterra integro- differential equations. Numer. Alogorithms 81, 57–77 (2023) Gander, W., Gautschi, W.: Adaptive quadrature revisited. BIT Numer. Math. 40, 84–101 (2000) Garrappa, R.: Numerical solution of fractional differential equations: a survey and a software tutorial. Mathematics 6(2), 16 (2018) Gottlieb, S., Gottlieb, D.: Spectral methods. Comput. Methods Appl. Mech. Eng. 4(9), 7504 (2009) Hans-Gorg, R., Martin, S., Lutz, T.: Robust numerical methods for singularly perturbed differential equa- tions: convection-diffusion-reaction and flow problems. Springer Ser. Comput. Math. 257, 1–604 (1994) Huseynov, I.T., Ahmadova, A., Fernandez, A., Mahmudov, N.I.: Explicit analytical solutions of incommen- surate fractional differential equation systems. Appl. Math. Comput. 390, 125590 (2021) Kadalbajoo,M.K., Patidar, K.C.: A survey of numerical techniques for solving singularly perturbed ordinary differential equations. Appl. Math. Comput. 130, 457–510 (2002) Li, C., Zeng, F., Liu, F.: Spectral approximation to the fractional integral and derivative. Fract. Calc. Appl. Anal. 15(3), 383–406 (2012) Logan, J.D.: An Introduction to Nonlinear Partial Differential Equations. Wiley, New York (1994) Motsa, S.: A new piecewise-quasilinearizationmethod for solving chaotic systems of initial value problems. Cent. Eur. J. Phys. 10(4), 936–946 (2012) Motsa, S.S., Dlamini, P., Khumalo, M.: A new multistage spectral relaxation method for solving chaotic initial value systems. Nonlinear Dyn. 72, 265–283 (2013) Ninghu, S.: Fractional Calculus for Hydrology, Soil Science and Geomechanics. Taylor & Francis, New York (2020) Nnakenyi, C.A., Weideman, J.A.C., Hale, N.: Spectral methods for fractional differential equations. AIMS Master Thesis (2015) Oloniiju, S.D., Goqo, S.P., Sibanda, P.: A Chebyshev spectral method for heat and mass transfer in MHD nanofluid flow with space fractional constitutive model. Front. Heat and Mass Transf. 13, 3–14 (2019) Oloniiju, S.D., Goqo, S.P., Sibanda, P.: A pseudo-spectral method for time distributed order two-sided space fractional differential equations. Taiwan. J. Math. 25(5), 959–979 (2021) Oloniiju, S.D., Mukwevho, N., Tijani, Y.O., Otegbeye, O.: Chebyshev pseudospectral method for fractional differential equations in non-overlapping partitioned domains. Appl. Math. 4(3), 950–974 (2024) Otegbeye, O.: On paired decoupled quasi-linearizationmethods for solving nonlinear systems of differential equations that model boundary layer fluid flow problems. Ph.D. thesis (2018) Owolabi, K., Atangana, A.: Numerical Methods for Fractional Differentiation. Springer, Berlin (2019) Peter, D., Martin, W.: Adaptive Numerical Solution of PDEs. De Gruyter, Berlin (2012) Podlubny, I.: Fractional Differential Equations. B Academic Press, San Diego. (1999) Qureshi, S., Akanbi, M.A., Shaikh, A.A., Wusu, A.A., Ogunlaran, O.M., Mahmoud, W., Osman, M.S.: A new adaptive nonlinear numerical method for singular and stiff differential problems. Alex. Eng. J. 74, 585–597 (2023) Ramos, J.I.: Linearization techniques for singularly-perturbed initial-value problems of ordinary differential equations. Appl. Math. Comput. 163, 1143–1163 (2005) Rossikhin, Y.A., Shitikova, M.V.: Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids. Appl. Mech. Rev. 50(1), 15–67 (1997) Samuel, F.M., Motsa, S.S.: Solving hyperbolic partial differential equations using a highly accurate mul- tidomain bivariate spectral collocation method. Wave Motion 88, 57–72 (2019) Sharma, K.K., Rai, R., Patidar, K.C.: A review on singularly perturbed differential equation with turning points and interior layers. Appl. Math. Comput. 219(3), 10575–10609 (2013) Shen, J., Tang, T.,Wang, L.-L.: SpectralMethods: Algorithms, Analysis andApplications, vol. 41. Springer, Berlin (2011) Sikora, B.: Remarks on the caputo fractional derivative. Matematyka Informatyka Na Uczelniach Tech- nicznych (MINUT) 5, 76–84 (2023) Sweilam,N.H.,Khader,M.M.: Efficient chebyshev spectralmethods for solvingmulti-term fractional orders differential equations. ANZIAM J. 51(4), 464–475 (2010) Trefethen, L.N.: Spectral Methods in MATLAB. Society for Industrial and Applied Mathematics, Provi- dence (2000) 123 100 Page 22 of 22 Journal of Nonlinear Science (2025) 35:100 Xu, Q., Hesthaven, J.S.: Stable multi-domain spectral penalty methods for fractional partial differential equations. J. Comput. Phys. 257, 241–258 (2014) Zayernouri, M., Em Karniadakis, G.: Fractional spectral collocation method. SIAM J. Sci. Comput. 36(1), A40–A62 (2014) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123 Adaptive Multidomain Numerical Solution for Singularly Perturbed Fractional Differential Equation: Chebyshev Pseudospectral Method Abstract 1 Introduction 2 Preliminaries and Definitions 3 Numerical Discretization via Pseudospectral Approximation 3.1 Approximation of the Function u(t) 3.2 Approximation in Equal Nonoverlapping Partitioned Domains 3.3 Adaptive Nonoverlapping Partitioned Domain for Initial-Value Problem 3.4 Error Estimates 4 Numerical Experimentation and Results 5 Conclusion References