In this study, the solution of the Hamilton-Jacobi equation (HJE) with holonomic Hamiltonian is investigated in terms of the first integrals of the corresponding Hamiltonian system. Holonomic functions are related to a specific type of partial differential equations called Pfaffian systems, whose solution space can be regarded as a finite-dimensional real vector space. In the finite-dimensional solution space, the existence of first integrals that define a solution of the HJE is characterized by a finite number of algebraic equations for finite-dimensional vectors, which can be easily solved and verified. The derived characterization was illustrated through a numerical example.