Accurate source localization of electroencephalographic (EEG) signals requires detailed information about the geometry and physical properties of head tissues. Indeed, these strongly influence the propagation of neural activity from the brain to the sensors. Finite difference methods (FDMs) are head modelling approaches relying on volumetric data information, which can be directly obtained using magnetic resonance (MR) imaging. The specific goal of this study is to develop a computationally efficient FDM solution that can flexibly integrate voxel-wise conductivity and anisotropy information. Given the high computational complexity of FDMs, we pay particular attention to attain a very low numerical error, as evaluated using exact analytical solutions for spherical volume conductor models. We then demonstrate the computational efficiency of our FDM numerical solver, by comparing it with alternative solutions. Finally, we apply the developed head modelling tool to high-resolution MR images from a real experimental subject, to demonstrate the potential added value of incorporating detailed voxel-wise conductivity and anisotropy information. Our results clearly show that the developed FDM can contribute to a more precise head modelling, and therefore to a more reliable use of EEG as a brain imaging tool.