In the present scientific article, an efficient operational matrix based on the famous Bernstein polynomials is applied for the numerical solution of two-dimensional non-linear variable order reaction-diffusion equation in porous media with given initial and boundary conditions. An operational matrix is constructed for fractional variable order differentiation w.r.to space variable x,y and time t, so that our proposed model is converted into a system of non-linear algebraic equations with the help of collocation method, which can be solved employing the Newton-Iteration method. The salient features of the article are finding the stability analysis and error bounds of the proposed method and also the validation and the effectiveness of the method through the RMS, L∞ and L2 errors. The physical presentation of the these errors for considered twodimensional non-linear variable order reaction-diffusion with their exact solutions shows the method is too good for finding the solution of these kind of problems.