We can not list the applications or the fields that use the anomalous sub-diffusion equations due to their wide area, one of these important applications are in the chemical reactions when a single substance tends to move from an area of high concentration to an area of low concentration until the concentration is equal across space. The mathematical model that describes these physical-chemical phenomena is called the reaction subdiffusion equation. In our study, we try to solve the 2D variable order version of these equations (2DVORSE) (linear and nonlinear) using an accurate numerical technique which is the weighted average finite difference method (WAFDM). We will study the stability of the resulting scheme using the fractional version of the John von Neumann stability analysis procedure. An accurate specific stability condition that is valid for some parameters in the resulting schemes is derived and checked. At the end of the study, we present some numerical examples to demonstrate the accuracy of the proposed technique.