<abstract><p>In this research article, we develop a powerful algorithm for numerical solutions to variable-order partial differential equations (PDEs). For the said method, we utilize properties of shifted Legendre polynomials to establish some operational matrices of variable-order differentiation and integration. With the help of the aforementioned operational matrices, we reduce the considered problem to a matrix type equation (equations). The resultant matrix equation is then solved by using computational software like Matlab to get the required numerical solution. Here it should be kept in mind that the proposed algorithm omits discretization and collocation which save much of time and memory. Further the numerical scheme based on operational matrices is one of the important procedure of spectral methods. The mentioned scheme is increasingly used for numerical analysis of various problems of differential as well as integral equations in previous many years. Pertinent examples are given to demonstrate the validity and efficiency of the method. Also some error analysis and comparison with traditional Haar wavelet collocations (HWCs) method is also provided to check the accuracy of the proposed scheme.</p></abstract>