We provide a novel discretization method with the help of the new operational matrices and fractal-fractional derivative operator for solving time fractal-fractional partial differential equations. To achieve our target, we consider the Bessel functions of the first kind to get the approximate solution with high precision. For the proposed problem, the basis functions together with their corresponding operational matrices are reduced to a system of algebraic equations. Besides, the error analysis of the method is thoroughly discussed. At last, to confirm the applicability and efficiency of the methodology, we implement several numerical tests. Furthermore, we discuss numerically HIV infection of the model of $\text{CD}{4}^{+}\text{T}$ cells.