Abstract-We propose and implement an algorithm based on reachability analysis to estimate the region of attraction (ROA) of an equilibrium point for nonlinear systems. The stability region is obtained via the computation of forward reachable sets. We compare our results with well-established techniques in this area. In particular, we consider the optimization of the Lyapunov function (LF) sub-level set using sum-of-squares (SOS) decomposition, and the computation of backward reachable sets of a target set using the viscosity solution of a time-dependant Hamilton-Jacobi-Isaacs (HJI) formulation. Our method can overcome many limitations imposed on the applicability of Lyapunov-based approaches, such as conservatism in estimating the stability region, and difficulties associated with choosing a suitable LF. This is due to the fact that our reachability algorithm does not require a LF in order to provide an estimate of the ROA. Various numerical examples show that our proposed approach can estimate the exact ROA quite accurately, and more importantly, scales moderately with the system dimension compared to alternative techniques.