In the present article, an efficient operational matrix based on the famous Laguerre polynomials is applied for the numerical solution of two-dimensional non-linear time fractional order reaction–diffusion equation. An operational matrix is constructed for fractional order differentiation and this operational matrix converts our proposed model into a system of non-linear algebraic equations through collocation which can be solved by using the Newton Iteration method. Assuming the surface layers are thermodynamically variant under some specified conditions, many insights and properties are deduced e.g., nonlocal diffusion equations and mass conservation of the binary species which are relevant to many engineering and physical problems. The salient features of present manuscript are finding the convergence analysis of the proposed scheme and also the validation and the exhibitions of effectiveness of the method using the order of convergence through the error analysis between the numerical solutions applying the proposed method and the analytical results for two existing problems. The prominent feature of the present article is the graphical presentations of the effect of reaction term on the behavior of solute profile of the considered model for different particular cases.