This paper presents an approximate method of solving a fractional (in the time variable) equation which describes the processes lying between heat and wave behavior. In most of applications related to fractional differential or fractional integro-differential equations, the numerical methods are limited to

$1+1$ (time + space) dimensions. This paper presents a different method for solving a fractional integro-differential equation which can handle more dimensional cases within a good approximation. The method is limited to cases when the initial condition is smooth enough with respect to space variables

$x$. In such cases the approach works also for

$(1+2)$ and

$(1+3$) dimensions. More precisely the approximation consists in the application of a finite subspace of an infinite basis in the time variable (Galerkin method) and discretization in space variables. In the final step, a large-scale system of linear equations with a non-symmetric matrix is solved with the use of the iterative GMRES method.