The gradient estimates in the generalized Orlicz space for weak solutions of a class of quasi-linear elliptic boundary value problems are obtained using the modern technique of extrapolation. The coefficients are assumed to have small BMO seminorms, and the boundary of the domain is sufficiently flat in the sense of Reifenberg. As a corollary, we apply our results to the variable Lebesgue spaces.