The problem of computing bounds on the region-of-attraction for systems with polynomial vector fields is considered. Invariant subsets of the region-of-attraction are characterized as sublevel sets of Lyapunov functions. Finite-dimensional polynomial parametrizations for Lyapunov functions are used. A methodology utilizing information from simulations to generate Lyapunov function candidates satisfying necessary conditions for bilinear constraints is proposed. The suitability of Lyapunov function candidates is assessed solving linear sum-of-squares optimization problems. Qualified candidates are used to compute invariant subsets of the region-of-attraction and to initialize various bilinear search strategies for further optimization. We illustrate the method on small examples from the literature and several control oriented systems. (c) 2008 Elsevier Ltd. All rights reserved,