This dissertation presents a decisive advance in the numerical solution and analysis of fractional diffusion, a relatively new but rapidly growing area of research. We exploit the cylindrical extension proposed and investigated by X. Cabre and J. Tan, in turn inspired by L. Caffarelli and L. Silvestre, to replace the intricate integral formulation of the fractional Laplacian, in a bounded domain, by a local elliptic PDE in one higher dimension with variable coefficient. Inspired in the aforementioned localization results, we propose a simple strategy to study discretization and solution techniques for problems involving fractional powers of elliptic operators. We develop a complete and rigorous a priori and interpolation error analyses. We also design and study an efficient solver, and develop a suitable a posteriori error analysis. We conclude showing the flexibility of our approach by analyzing a fractional space-time parabolic equation.