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 Vn in a rational functional form approximating a maximal Lyapunov function VM 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 \textit{Mathematica}-package to find such Lyapunov functions Vn for a class of hybrid (piecewise non-linear) systems, where the dynamics is continuous on the boundary of the different regimes in the state space. In addition, a computationally feasible method is proposed to estimate the DOA using a maximal fitting hypersphere.