Abstract-We explore region of attraction (ROA) estimation for polynomial systems via the numerical solution of polynomial equations. Computing an optimal, stable sub-level set of a Lyapunov function is first posed as a polynomial optimization problem. Solutions to this optimization problem are found by solving a polynomial system of equations using techniques from numerical algebraic geometry. This system describes KKT points and singular points not satisfying a regularity condition. Though this system has exponentially many solutions, the proposed method trivially parallelizes and is practical for problems of moderate dimension and degree. In suitably generic settings, the method solves the underlying optimization problem to arbitrary precision, which could make it a useful tool for studying popular semidefinite programming based relaxations used in ROA analysis.