An implementation of the localized boundary-domain integral-equation (LBDIE) method for the numerical solution of the Neumann boundary-value problem for a second-order linear elliptic PDE with variable coefficient is discussed. The LBDIE method uses a specially constructed localized parametrix (Levi function) to reduce the BVP to a LBDIE. After employing a mesh-based discretization, the integral equation is reduced to a sparse system of linear algebraic equations that is solved numerically. Since the Neumann BVP is not unconditionally and uniquely solvable, neither is the LBDIE. Numerical implementation of the finite-dimensional perturbation approach that reduces the integral equation to an unconditionally and uniquely solvable equation, is also discussed
