Lyapunov functions are commonly involved in the analysis of the stability of linear and nonlinear dynamical systems. Despite the fact that there is no generic procedure for creating these functions, many authors use polynomials in p-forms as candidates for constructing Lyapunov functions, while others restrict the construction to quadratic forms. We proposed a method for constructing polynomial Lyapunov functions that are not necessary in a form by focusing on the positive and negative definiteness of the Lyapunov candidate and the Hessian of its derivative, as well as employing the sum of square decomposition. The idea of Newton polytopes was used to transform the problem into a system of algebraic equations that were solved using the polynomial homotopy continuation method. Our method can produce several possibilities of Lyapunov functions for a given candidate. The sample test conducted demonstrates that the method developed is promising.