In this work, we study, from both analytical and numerical points of view, a heat conduction model which is based on the Moore-Gibson-Thompson equation. The second gradient effects are also included. First, the existence of a unique solution is proved by using the theory of linear semigroups, and the exponential energy decay is also shown when the constitutive tensors are homogeneous. The analyticity of the semigroup is also discussed in the isotropic case, and its spatial behavior is studied. The spatial exponential decay is also proved. Then, we provide the numerical analysis of a fully discrete approximation obtained by using the finite element method and an implicit Euler scheme. A discrete stability property is shown, and some a priori error estimates are derived, from which the linear convergence is concluded under suitable regularity conditions. Finally, some one-dimensional numerical simulations are presented to demonstrate the accuracy of the approximations and the behavior of the discrete energy.