In this paper, a compact finite difference scheme with global convergence order O(τ+h) is derived for fractional sub-diffusion equations with the spatially variable coefficient subject to Neumann boundary conditions. The difficulty caused by the variable coefficient and the Neumann boundary conditions is overcome by subtle decomposition of the coefficient matrices. The stability and convergence of the proposed scheme are studied using its matrix form by the energy method. The theoretical results are supported by numerical experiments.