In recent years, random functional or stochastic equations have been reported in a large class of problems. In many cases, an exact analytical solution of such equations is not available and, therefore, is of great importance to obtain their numerical approximation. This study presents a numerical technique based on Bernstein operational matrices for a family of stochastic fractional integro-differential equations (SFIDE) by means of the trapezoidal rule. A rele- vant feature of this method is the conversion of the SFIDE into a linear system of algebraic equations that can be analyzed by numerical methods. An upper error bound, the convergence, and error analysis of the scheme are investigated. Three examples illustrate the accuracy and performance of the technique