In this study, an effective numerical approach based on the hybrid of block-pulse and parabolic functions (PBPFs) is suggested to obtain an approximate solution of a system of nonlinear stochastic Itˆo-Volterra integral equations of fractional order. For this aim, we first introduce these functions and express some of their properties and then calculate fractional and stochastic operational matrices of integration based on these functions. Using the properties of PBPFs and obtained operational matrices, the system of nonlinear stochastic Itˆo-Volterra integral equations of fractional order convert to a nonlinear system of algebraic equations which can be easily solved by using Newton’s method. Moreover, in order to show the rate of convergence of the suggested approach, we present several theorems on convergence analysis and error estimation which demonstrate the rate of convergence of the proposed method for solving this nonlinear system is O(h3). Finally, two examples are included to illustrate the validity, applicability and efficiency of the proposed technique.