Edge turbulence models will be crucial to interprete data and help design future nuclear fusion experiments and power plants, like ITER and EU-DEMO. There is a need to develop numerical methods that can effectively scale with the reactor dimensions and handle the highly anisotropic nature of turbulence. As an initial step, we analyze a system describing the dynamics of a fast wave commonly observed in plasma physics: the electrostatic shear Alfvén waves (SAWs). We propose a finite difference method on staggered grids for wave-like problems, which is rigorously energy-preserving and accurate. This method emulates certain aspects of the summation-by-parts (SBP) operator framework, particularly the preservation of the divergence theorem at the discrete level. The design of this method is intended to be versatile, making it suitable for a wide range of wave problems characterized by divergence-free velocity fields. As a subsequent step, we incorporate the ExB-drift dynamics, which evolves on a characteristic time scale slower than that of the waves. In this system, SAWs represent the fastest phenomena parallel to the magnetic field, while the E×B-drift dynamics predominantly occurs perpendicular to the equilibrium magnetic field. To exploit this time-scale separation, we integrate the SAWs using an unconditionally stable implicit method, while the ExB-drift dynamics is treated explicitly. Temporal integration is achieved using the globally stiffly accurate IMEX Runge-Kutta, called BPR((3,5,3)) scheme. Consequently, the overall time step of the algorithm is determined by the slower E×B-dynamics. By combining stable spatial discretization for SAWs and an IMEX approach for time integration, we achieve both stability and computational efficiency in our simulations.