Helmholtz equation is a class of elliptic partial differential equations that describe electromagnetic waves, and are widely used in mechanics, acoustics, electromagnetism, and other fields. In order to eliminate the pollution effect of high wavenumbers, the traditional method for numerically solving Helmholtz equation is to refine the grid, which not only increases the time complexity, but also usually makes the discrete matrix ill conditioned. Therefore, it is necessary to find an efficient numerical method for any wavenumbers. Based on the finite volume method, variable limit factors are introduced to completely convert the differential equations into integral equations. A discrete scheme containing a tridiagonal matrix is constructed using the univariable three-point and bivariable nine-point Lagrange interpolation formulas to perform numerical solutions of one-dimensional and two-dimensional Helmholtz equations using the variable limit integration method, respectively. The proposed method is suitable for arbitrary wave numbers, and the physical meaning of the solution process is clear. For the one-dimensional Helmholtz equation, the influence of the variable limit factor on the error is studied. The error estimation of the numerical solution is performed using Taylor expansion and Lagrange interpolation residual formula, and it is proved that the truncation error of the discrete scheme reaches second order. Numerical examples show that when the variable limit factor and step size of the discrete scheme are equal, the error order is lower. For the two-dimensional Helmholtz equation, the influence of different wave numbers on the numerical solution is investigated. It is proved that the truncation error of the discrete scheme reaches third order. Numerical examples indicate that the numerical scheme has good accuracy for different wave numbers, and high wavenumbers do not cause the pollution effect.