Objective: This paper proposes a stable volume integral equation (VIE) solver based on polarization/magnetization currents, for the accurate and efficient computation of the electromagnetic scattering from highly inhomogeneous and high contrast objects. Methods: We employ the Galerkin Method of Moments to discretize the formulation with discontinuous piecewise linear basis functions on uniform voxelized grids, allowing for the acceleration of the associated matrix-vector products in an iterative solver, with the help of FFT. Results: Numerical experiments are conducted to study the accuracy and convergence properties of the proposed framework. Their results are compared against standard low order (piecewise constant) discretization schemes, a more conventional VIE formulation based on electric flux densities, and a commercial software package that employs the finite difference time domain method. Conclusion: The results illustrate the superior accuracy and wellconditioned properties of the proposed scheme. Significance: The developed solver can be applied to accurately analyze complex geometries, including realistic human body models, typically used in modeling the interactions between electromagnetic waves and biological tissue, that arise in magnetic resonance scanners.