      The topic of fractional differential equations (FDEs) has become an active area of research over the past several decades. The study of existence and uniqueness of solutions, and of the evolution of systems described by FDEs is of theoretical, as well as practical interest to mathematicians and scientists who aim to model the behaviour of complex dynamical systems. The main advantage of fractional calculus operators is in their ability to capture non-local and long-term memory effects [1]}, [2]}. This property allows the development of more realistic models using FDEs for complex phenomena, such as anomalous diffusion, the behaviour of viscoelastic materials, transport properties, and fluid flows [3]}. Generally, real-world processes are non-linear, and thus described by FDEs containing non-linearities, the exact solutions to which are often not possible to obtain. This has motivated the development of approximate solution methods, such as the numerical-analytic method, which combines deriving an approximate solution in analytic form with the numerical calculation of the parameters describing the solution's behaviour.
