Let $\alpha \in (0, \infty)$ satisfy $n-1 < \alpha < n$ for some $n \in \mathbb{N}$. The Caputo fractional derivative of order $\alpha$ for a function $f \in \mathcal{C}^n(0, \infty)$ is defined as
$$
D_t^{\alpha} f(t) = \frac{1}{\Gamma(n - \alpha)} \int_0^t \frac{f^{(n)}(\tau)}{(t - \tau)^{1 - n + \alpha}} \, d\tau.
$$
Similarly, for $\alpha \in (0, \infty)$ and $f \in \mathcal{C}^0(0, \infty)$, the Riemann--Liouville fractional integral of order $\alpha$ is given by
$$
I_t^{\alpha} f(t) = \frac{1}{\Gamma(\alpha)} \int_0^t \frac{f(\tau)}{(t - \tau)^{1-\alpha}} \, d\tau.
$$
In this talk, we present a pseudospectral approach to construct operational matrices based on shifted Chebyshev polynomials for the numerical approximation of Caputo fractional derivatives and Riemann--Liouville fractional integrals. To ensure numerical stability in the construction of these matrices, we employ variable precision arithmetic. We then apply the resulting Caputo differentiation matrices to solve Caputo-type advection--diffusion equations in one and multiple spatial dimensions, where the discretization reduces to a Sylvester (tensor) equation. Numerical experiments involving highly oscillatory time-dependent functions demonstrate the accuracy, stability, and effectiveness of the proposed method.