The superior ability of plasmonic structures to manipulate light has propelled their extensive applications in nanophotonics techniques and devices. Computational electromagnetics plays a critical role in characterizing and optimizing the nanometallic structures. In this paper, a general numerical algorithm, which is different from the commonly used discrete dipole approximation, the finite-difference time-domain, and the surface integral equation (SIE) method, is proposed to model plasmonic nanostructures. In this algorithm, the generalized impedance boundary condition (GIBC) based on the finite element method (FEM) is formulated and converted to the SIE. The plasmonic nanostructures with arbitrary inhomogeneity and shapes are modeled by the FEM. Their complex electromagnetic interactions are accurately described by the SIE method. As a result, the near field of plasmonic nanostructures can be accurately calculated. The higher order basis functions, together with the multifrontal massively parallel sparse direct solver, are involved to provide a higher order accurate and fast solver. © 2011 IEEE.