A numerical method is proposed for modeling electron cyclotron resonance heating (ECRH) and current drive in fusion devices with, generally speaking, nonlinear wave–particle interaction. Considered is the case of second harmonic extraordinary wave resonance. The method is applied to model the electron distribution function in a simplified magnetic field geometry. Various combined effects of wave–particle interaction and Coulomb collisions are demonstrated. The model is compared to the bounce-averaged quasilinear Fokker–Planck equation for the case of perpendicular propagation of the wave beam with respect to the magnetic field. In addition, the absorption coefficients computed from these two models are compared to the results of both the linear and the adiabatic model. Significant differences between results are found for present-day ECRH power levels and fusion device parameters.