An unsteady problem is considered for a space-fractional diffusion equation in a bounded domain. A first-order evolutionary equation containing a fractional power of an elliptic operator of second order is studied for general boundary conditions of Robin type. Finite element approximation in space is employed. To construct approximation in time, regularized two-level schemes are used. The numerical implementation is based on solving the equation with the fractional power of the elliptic operator using an auxiliary Cauchy problem for a pseudo-parabolic equation. The results of numerical experiments are presented for a model two-dimensional problem.