This paper presents an efficient polynomial kernel-based method for solving multi-term nonlinear time-fractional diffusion systems in high-dimensional regular and irregular domains. Time-fractional diffusion equations are increasingly recognized for their ability to model complex phenomena across various disciplines, including physics, engineering, and finance. By incorporating fractional-order temporal derivatives, these equations effectively capture nonlocal memory effects and intricate dynamics in the studied systems. We aim to provide an effective method for obtaining suitable approximate solutions to the time-fractional diffusion system of equations, which presents a significant challenge in numerical methods. We utilized the properties of Bernoulli polynomials to construct approximations in a finite-dimensional reproducing kernel space. The time-fractional system of equations is discretized using a polynomial kernel-based technique in the spatial direction. Then, an accurate high-order backward differentiation formula is utilized for time discretization. A thorough convergence analysis is conducted, establishing error bounds for the proposed method. To validate the versatility and effectiveness of the method, we conduct several numerical simulations demonstrating its performance across two and three dimensional domains.
