Electromagnetic metamaterials are artificial media assembled from components which have dimensions much smaller than the wavelength of the illuminating radiation. It has been demonstrated that these media can have properties beyond those found in conventional materials with important applications to many areas of science and engineering. Among the collection of such materials currently garnering significant attention are the Hyperbolic Metamaterials. These are highly anisotropic structures which have a hyperbolic dispersion relation due to the fact that one principal component of the relative permittivity or permeability tensor has the opposite sign of the other two. For such materials scattered wave information at all scales is transmitted far away which suggests a number of important applications, the most immediate of which is imaging objects below the diffraction limit (superlensing). Clearly it is of great current interest to have numerical simulation capabilities for these fascinating materials. As the hyperbolic response is quite strong and its nature is very sensitive, numerical simulations of these configurations should be robust and highly accurate. For this reason we focus on High–Order Spectral algorithms which efficiently produce high fidelity solutions. More specifically, we describe a High–Order Perturbation of Surfaces approach which enjoys the greatly reduced operation counts and memory savings of interfacial methods while avoiding the complexities and indefinite linear systems faced by Integral Equation algorithms. We give a full discussion of our formulation in terms of Impedance–Impedance Operators (which avoids spurious singularities of other formulations) and implementation details, followed by numerical validation and simulation of results which appear in the engineering literature.