This study is concerned with numerical solutions of a Magneto-Thermo-Elasticity (MTE) model via a combination of energy-conserving finite difference method (FDM) with Richardson extrapolation methods (REMs). Firstly, by introducing two auxiliary functions and using second-order centered FDM and Crank-Nicolson method to approximate spatial and temporal derivatives, respectively, a two-level energy-conserving FDM is established for a MTE model. The priori estimation, solvability, and convergence are derived rigorously by using the discrete energy method. Secondly, to improve computational efficiency, a class of REMs are also designed by constructing the symbolic expansion of numerical solutions. Finally, numerical results confirm the efficiency of the proposed algorithms and the exactness of the theoretical findings.