An integral equation (IE) based field solver to compute the scattered fields from spatially dispersive metasurfaces is proposed and numerically confirmed using various examples involving physical unit cells. The work is a continuation of Part-1 [1], which proposed the basic methodology of representing spatially dispersive metasurface structure in the spatial frequency domain, k||. By representing the angular dependence of the surface susceptibilities in k|| as a ratio of two polynomials, the standard generalized sheet transition conditions (GSTCs) have been extended to include the spatial derivatives of both the difference and average fields around the metasurface. These extended boundary conditions are successfully integrated here into a standard IE-GSTC solver, which leads to the new IE-GSTC-SD simulation framework. The proposed IE-GSTC-SD platform is applied to various uniform metasurfaces, including a practical short conducting wire unit cell, as a representative practical example, for various cases of finite-sized flat and curvilinear surfaces. In all cases, computed field distributions are successfully validated, either against the semi-analytical Fourier decomposition method or the brute-force full-wave simulation of volumetric metasurfaces in the commercial Ansys FEM-HFSS simulator.