<p>This study focuses on developing a numerical technique to approximate the solution of a variable-order time-fractional (VO TF) integro-partial differential equation (IPDE) with a weakly singular kernel. The numerical technique proposed in this paper is based on two-dimensional (2D) wavelets, namely fractional-order Chelyshkov wavelets (FO CWs) and Bernoulli wavelets (BWs). The main novelty lies in expressing any square-integrable function in terms of FO CWs and BWs, along with the corresponding integral and variable-order operational matrices (OMs), which are constructed using these wavelets to facilitate and enhance the proposed numerical approach. With the help of these matrices and suitable collocation points, the given problem is reduced to an algebraic system of equations (ASEs), which is then solved using Newton’s iteration method. Additionally, a variable-order time-fractional integro-partial differential equation with a weakly singular kernel is solved using the spectral collocation method (SCM) based on SL–GL–C and SC–GR–C. A comparison of the numerical results shows that the two-dimensional wavelet method provides significantly more accurate results than the SCM. Therefore, the proposed two-dimensional methods are more effective than the SCM for solving variable-order time-fractional integro-partial differential equations. In addition, several relevant theorems are discussed, including the analysis of convergence and error estimates for the 2D wavelet method. Several examples are presented to illustrate the effectiveness and applicability of the proposed method.</p>