The main purpose of this book is to study the following problem: If a system of partial differential equation decays polynomially to zero, what are the conditions for which its numerical approximation still decreases polynomially to zero uniformly with respect to the step size? To answer this question, we consider two different methods: The first one called the resolvent method is to link the decay of the system to equilibrium to the asymptotic behavior of the resolvent along the imaginary axis. The second one called the energy method or the renorming method consists in choosing an appropriate Lyapunov functional, to prove that the total energy of solutions associated to the family of the studied systems is decreasing with a decay rate given explicitly. In order to illustrate these results, we consider examples from the theory of thermoelasticity: a thermoelastic wave system and a thermoelastic beam system. For a better understanding of this asymptotic behavior numerical simulations are also presented.