In this paper we analyze local asymptotic stability of switched hybrid systems, whose subsystems have polynomial vector fields, by discovering multiple Lyapunov functions in quadratic forms. We start with an algebraizable sufficient condition for the existence of quadratic multiple Lyapunov functions. Then, since different discrete modes are considered, we apply real root classification together with a projection operator to underapproximate this sufficient condition step by step, arriving at a set of semialgebraic sets which only involve the coefficients of the preassumed multiple Lyapunov function. Note that for each step, we additionally use the information on the different discrete modes to optimize our intermediate computation results. Afterward, we compute a sample point in the corresponding semialgebraic set for the coefficients, resulting in a multiple Lyapunov function. Finally, we test our approach on some examples using a prototypical implementation and compare it with the generic quantifier elimination based method and the sum of squares based method. These computation and comparison results show the applicability and promise of our approach. Furthermore, our approach is extendable for piecewise affine systems and nonpolynomial switched systems by discovering multiple Lyapunov functions beyond quadratic forms.