Based on a recent duality theory for linear differential inclusions (LDIs), the condition for stability of an LDI in terms of one Lyapunov function can be easily derived from that in terms of its conjugate function. This paper uses a particular pair of conjugate functions, the convex hull of quadratics and the maximum of quadratics, for the purpose of estimating the domain of attraction for systems with saturation nonlinearities. To this end, the nonlinear system is locally transformed into a pararnetertized LDI system with an effective approach which enables optimization on the parameter of the LDI along with the optimization of the Lyapunov functions. The optimization problems are derived for both the convex hull and the max functions, and the domain of attraction is estimated with both the convex hull of ellipsoids and the intersection of ellipsoids. A numerical example demonstrates the effectiveness of this paper's methods. (c) 2005 Elsevier Ltd. All rights reserved.