Modeling the propagation and damping of electromagnetic waves in a hot magnetized plasma is difficult due to spatial dispersion. In such media, the dielectric response becomes non-local and the wave equation an integro-differential equation. In the application of RF heating and current drive in tokamak plasmas, the finite Larmor radius (FLR) causes spatial dispersion, which gives rise to physical phenomena such as higher harmonic ion cyclotron damping and mode conversion to electrostatic waves. In this paper, a new numerical method based on an iterative wavelet finite element scheme is presented, which is suitable for adding non-local effects to the wave equation by iterations. To verify the method, we apply it to a case of one-dimensional fast wave heating at the second harmonic ion cyclotron resonance, and study mode conversion to ion Bernstein waves (IBW) in a toroidal plasma. Comparison with a local (truncated FLR) model showed good agreement in general. The observed difference is in the damping of the IBW, where the proposed method predicts stronger damping on the IBW.