This article tackles the problem of estimating the domain of attraction of a Lur'e system, that is the feedback interconnection of a linear time-invariant system with a memoryless static operator. When the dimension of the system is large, numerical approaches based on simulations become prohibitive from a computational point of view. On the other hand, classical analytical techniques based on Lyapunov functions provide conservative estimates because they usually consider quadratic Lyapunov functions. Another limit is given by the fact that they deal with contractively invariant sets, which are sets where the derivative of the Lyapunov function along the trajectories is negative. Methods to reduce their conservativeness are still a challenging subject of research and desirable for many practical applications. In this paper, we try to combine the information given by more Lyapunov functions together in order to enlarge the estimate of the domain of attraction. The novelty of our approach lies in the fact that the sets we are considering are invariant but not necessarily contractively invariant. We assume that an estimate of the attraction domain is already known. In order to show that a set is part of the attraction domain, it is sufficient to prove that all the trajectories starting from it reach the current estimate in a finite time. This concept provides a method to iteratively improve the attraction domain estimate using different Lyapunov functions without limiting the analysis to contractively invariant sets. After developing a general theory, we only resort to the use of quadratic Lyapunov functions because of the computational appeal given by LMI solvers.