This article presents the construction and MATLAB implementation of the virtual element method for solving a time-fractional convection-diffusion equation characterized with fractional-order of \(\alpha \in (0,1)\) in temporal direction. The study of the following fractional problem over regular Voronoi meshes is carried out to demonstrate the implementation. \(\begin{aligned} ^{R}D^{\alpha }_{t}u(\varvec{x},t)-\Delta u (\varvec{x},t)+\varvec{b}\cdot \nabla {u}(\varvec{x},t)=f(\varvec{x},t)~\text {in}~\varvec{x}\in \Omega ,~t\in (0,T], \end{aligned}\) where \(\Omega\) represents a \(\mathbb {R}^2\) domain in space, \(\alpha\) the fractional-order, and t is temporal variable. The fractional Grünwald-Letnikov approximation serves as the foundation for our methodology. Using the energy projection operator and the discrete maximal regularity, we develop a discrete scheme that maintains both polynomial consistency and stability by definition. Complete details of implementation and original codes are included for solution of fully discrete scheme. The practical usefulness of the suggested method is highlighted by numerical results that validate the convergence in the \(L^2\) -norm and \(H^1\) -seminorm over regular Voronoi mesh configuration.
