This paper investigated a numerical method for a two-dimensional variable coefficient evolution equation (VCEE), utilizing the alternating direction implicit (ADI) method and an extrapolation formula. The time derivative was discretised by the backward Euler (BE) scheme on a uniform mesh and the finite difference method (FDM) was applied to spatial discretization. We proved an priori estimate and the error bound of the solution to the difference scheme using the energy analysis method, and verified the uniqueness, stability, and convergence of the proposed scheme. To further improve numerical accuracy, we introduced a Richardson extrapolation method, which enhances the global accuracy to fourth order. Finally, some numerical examples were provided to demonstrate the validity of the theoretical analysis.