We present a novel formulation, based on the latest advancement in polynomial system solving via linear algebra, for identifying limit cycles in general n-dimensional autonomous nonlinear polynomial systems. The condition for the existence of an algebraic limit cycle is first set up and cast into a Macaulay matrix format whereby polynomials are regarded as coefficient vectors of monomials. This results in a system of polynomial equations whose roots are solved through the null space of another Macaulay matrix. This two-level Macaulay matrix approach relies solely on linear algebra and eigenvalue computation with robust numerical implementation. Furthermore, a state immersion technique further enlarges the scope to cover also non-polynomial (including exponential and logarithmic) limit cycles. Application examples are given to demonstrate the efficacy of the proposed framework.