We propose a fast and stable method for constructing matrix approximations to fractional integral operators applied to series in Chebyshev fractional polynomials. Based on a recurrence relation satisfied by the definite integrals of mapped Chebyshev polynomials with a fractional weight, the proposed method significantly outperforms existing approaches. Through numerical examples, we highlight the broad applicability of these matrix approximations, including the solution of boundary value problems for fractional integral and differential equations. Additional applications include fractional differential equation initial value problems and fractional eigenvalue problems.