In this article, we study the numerical technique for variable‐order fractional reaction‐diffusion and subdiffusion equations that the fractional derivative is described in Caputo's sense. The discrete scheme is developed based on Lucas multiwavelet functions and also modified and pseudo‐operational matrices. Under suitable properties of these matrices, we present the computational algorithm with high accuracy for solving the proposed problems. Modified and pseudo‐operational matrices are employed to achieve the nonlinear algebraic equation corresponding to the proposed problems. In addition, the convergence of the approximate solution to the exact solution is proven by providing an upper bound of error estimate. Numerical experiments for both classes of problems are presented to confirm our theoretical analysis.