A new Drude two critical points (D-2CP) multiresolution time domain numerical method for the simulation of surface plasmon polariton in metallic inclusions is presented. A D-2CP model is used to accurately describe the dielectric function of the metallic inclusions, and it is compared with other Drude models reported in the literature. The superior accuracy given by the D-2CP model avoids the formation of spurious reflections at the boundary interface between metal and dielectric, which can greatly affect the accuracy of the numerical results.