This paper considers the problem of approximating the "maximal"region of attraction (the set that contains all asymptotically stable sets) of any given set of locally exponentially stable nonlinear Ordinary Differential Equations (ODEs) with a sufficiently smooth vector field. Given a locally exponential stable ODE with a differentiable vector field, we show that there exists a globally Lipschitz continuous converse Lyapunov function whose 1-sublevel set is equal to the maximal region of attraction of the ODE. We then propose a sequence of d-degree Sum-of-Squares (SOS) programming problems that yields a sequence of polynomials that converges to our proposed converse Lyapunov function uniformly from above in the L1 norm. We show that each member of the sequence of 1-sublevel sets of the polynomial solutions to our proposed sequence of SOS programming problems are certifiably contained inside the maximal region of attraction of the ODE, and moreover, we show that this sequence of sublevel sets converges to the maximal region of attraction of the ODE with respect to the volume metric. We provide numerical examples of estimations of the maximal region of attraction for the Van der Pol oscillator and a three dimensional servomechanism.