In this article a method is presented to find systematically the domain of attraction (DOA) of hybrid non-linear systems. It has already been shown that there exists a sequence of special kind of Lyapunov functions $V_n$ in a rational functional form approximating a maximal Lyapunov function $V_M$ that can be used to find an estimation for the DOA. Based on this idea, an improved method has been developed and implemented in a Mathematica-package to find such Lyapunov functions $V_n$ for a class of hybrid (piecewise...