In this paper, we study the finite volume discretization method for balanced fractional diffusion equations where the fractional differential operators are comprised of both Riemann-Liouville and Caputo fractional derivatives. The main advantage of this approach is that new symmetric positive definite Toeplitz-like linear systems can be constructed for solving balanced fractional diffusion equations when diffusion functions are non-constant. It is different from non-symmetric Toeplitz-like linear systems usually obtained by existing numerical methods for fractional diffusion equations. The preconditioned conjugate gradient method with circulant and banded preconditioners can be applied to solve the proposed symmetric positive definite Toeplitz-like linear systems. Numerical examples, for both of one- and two- dimensional cases, are given to demonstrate the good accuracy of the finite volume discretization method and the fast convergence of the preconditioned conjugate gradient method.