This study investigates the impact of mesh refinement and polynomial interpolation order on the accuracy of finite element method (FEM) simulations of temporal interference (TI) brain stimulation. A multilayer spherical head model was employed, and an analytical solution was derived to compute the relative error in the estimated maximum modulation amplitude within a predefined region of interest (ROI). Six mesh configurations were compared, combining uniform and locally refined meshes with linear (P1) and quadratic (P2) elements. Results show that local refinement near current injection sites significantly improves accuracy, even when using fewer total elements. Moreover, P2 elements consistently outperform P1, particularly in configurations with short electrode separations, where field gradients are more intense. Tools such as FEniCS facilitate the selection of the polynomial interpolation order, enabling straightforward comparison between different finite element spaces. Visualizations of modulation amplitude and absolute error across different configurations confirmed that the largest errors occur near singularities induced by point current sources and across tissue interfaces. Future work will extend this analysis to anatomically realistic head geometries and explore further comparisons between complete electrode models and point-source approximations using higher-order elements.