In this article, we derive a computational scheme joining the finite difference and operational matrix methods to numerically simulate fractional differential equations of distributed order with Riesz derivatives subject to the initial and Dirichlet boundary conditions. The finite difference method is used to discretize Riesz‐space fractional derivatives, and the operational matrix method is based on fractional shifted Gegenbauer functions for the approximation of time derivatives. The derived method transforms the given equation into a system of linear algebraic equations. For the numerical solution, we derive the optimal error bound and prove the theoretical unconditional stability with respect to the
‐norm. Furthermore, we numerically verify the stability of the proposed method. We observe the accuracy of the second order of the given schemes. The accuracy and effectiveness of the method are checked to solve some problems of telegraph equations. The CPU running time for the method with a fractional shifted Gegenbauer basis is much less than for methods with other bases.