We study the stability of a two-dimensional fractional reaction-diffusion system under the Caputo differential operator in time. Our model is based on the Grey-Scott model, a well-known coupled reaction-diffusion system that describes the interaction between two chemical species. The diffusion term captures the species special spread, while the nonlinear term in the system describes the chemical reaction, resulting in a wide range of difficult, self-organizing patterns, including spots, stripes, or spirals, depending on the parameter values. We derive conditions for local stability of the homogeneous equilibrium by linearizing the system and analyzing the eigenvalues of the Jacobian. Furthermore, we construct appropriate Lyapunov functionals to establish global asymptotic stability of the discrete model under suitable conditions. This approach seeks to provide a robust framework for analyzing complex dynamical behaviors in systems governed by fractional-order in-time reaction-diffusion systems. The numerical simulations employ the Chebyshev spectral method for spatial discretization and the $L1$ scheme for fractional time derivatives. These simulations validate the theoretical findings, demonstrating the model's ability to replicate intricate patterns often observed in reaction-diffusion systems. The results suggest that the fractional-order framework enhances the understanding of pattern formation in such systems, making this model a valuable tool for studying anomalous diffusion and non-local dynamics in biological and chemical processes.

Nonlinear stochastic Itô–Volterra integral equations (NSIVIEs) represent systems whose current state is influenced by random fluctuations and is dependent on previous information. These equations appear in many real-world scenarios, including engineering systems, biological processes, financial markets, heterogeneous media, complex transport phenomena, and viscoelastic materials. Strong numerical frameworks are required because analytical solutions for these equations are rarely available, particularly when nonlinearities and random fluctuations are present. To effectively solve NSIVIEs, in this study we propose a new hybrid numerical framework that combines Monte Carlo simulation and Legendre spectral collocation. By using orthogonal polynomial basis functions to approximate the solution, this method provides spectral accuracy while handling the hereditary memory component of the Volterra equation through a high-order Legendre spectral collocation method. A precise statistical treatment of the random fluctuations is made possible by simultaneously addressing the stochastic Itô noise through Monte Carlo sampling across numerous independent realizations. We perform a thorough convergence analysis and obtain explicit error bounds that measure the decrease in approximation error with increasing spectral resolution and Monte Carlo sample count. Numerical experiments show that the method can accurately reproduce complex stochastic behaviors and validate theoretical predictions.