## 1 Introduction  
Partial differential equations (PDEs) are widely used in physics and engineering, with applications including heat transfer, wave propagation, electromagnetism, and fluid dynamics. They are also applied in economics to calculate price alternatives and in biology to simulate population dynamics and disease propagation. In medicine, PDEs play an important role in imaging. Furthermore, they are used in geosciences to predict weather and simulate groundwater movement and in material science to study chemical reactions and diffusion. Their broad applications make them critical for comprehending and forecasting complex, dynamic systems in theoretical and practical settings. For some PDE applications, one can consult [1, 2]. Many numerical approaches were used to solve different types of PDEs. For example, in [3], the authors solved PDEs using deep learning and physical constraints. In [4], the authors presented some numerical algorithms for handling some non-linear PDEs. In [5], a variational quantum algorithm was followed for handling PDEs. Some approximations for certain stochastic parabolic PDEs were presented in [6]. In [7], a boundary element method was followed to solve one-dimensional nonlinear parabolic PDEs. In [8], the authors proposed a proper orthogonal decomposition method for handling certain PDEs. In [9], a numerical approach for a heat transfer model was proposed. In [10], the authors used a tau approach using Lucas polynomials for the time-fractional diffusion equation. Spectral methods are essential in numerical analysis due to their advantages over other numerical approaches (see, for example, [11; 12]). Approximate solutions derived from these methods are frequently represented as polynomial combinations, indicating their close relationship with special functions, especially orthogonal polynomials. There are three crucial spectral methods: collocation, tau, and the Galerkin method. All methods have their benefits and characteristics. The choice of a suitable spectral method is based on the type of differential equation studied and the underlying conditions (see, for example, [13; 14; 15; 16]). For other contributions concerning spectral methods, one can refer to [17; 18; 19]. Matrix representations of differentiation operations on a collection of basis functions are known as operational matrices of derivatives (OMDs) and are often employed in numerical approaches to solve various DEs. These matrices are effective tools to handle these DEs. The transformation of ordinary DEs into an algebraic system is possible by expressing the derivative operator as an operational matrix. It is useful in the collocation method ([20; 21]). Furthermore, when dealing with FDEs, the operational matrices of fractional derivatives will be the backbone for designing algorithms for these FDEs. For instance, the authors of [22] used shifted Legendre-Laguerre operational matrices to treat certain time-delay fractional diffusion equations. In [23], the author derived generalized Jacobi Galerkin OMDs and used them to handle specific FDEs. Other Legendre matrices were used in [24]. Operational matrices based on hat functions were used in [25] to treat some integral DEs. Galerkin OMDs were used in [26] to treat some multi-dimensional FDEs. In [27], OMDs of some finite classes of orthogonal polynomials were employed to handle FDEs. Other contributions regarding the use of different OMDs can be found in [28; 29; 30]. In both theoretical and applied mathematics, Chebyshev polynomials (CPs) are widely used; see [31; 32]. Approximation theory and numerical analysis are two areas where CPs are important. Due to their orthogonality, these polynomials provide other desirable characteristics, including minimization of error and oscillation (see, for example, [33; 34]). CPs are related to the classical Jacobi polynomials. Four kinds of CPs are special ones. All of these kinds have their characteristics and applications.