In the current paper, we introduce an efficient methodology to solve nonlinear stochastic differential equations (SDEs) driven by variable order fractional Brownian motion (vofBm) with appropriate initial condition. SDEs have many applications throughout pure mathematics and are used to model various behaviors of stochastic models such as stock prices, random growth models or physical systems that are subjected to thermal fluctuations. The mechanism of our proposed approach, which is based on Bernoulli polynomials operational matrices, is that it first transforms the problem under consideration into a nonlinear stochastic integral equation (SIE) driven by vofBm by using given initial condition and integrating from both sides of nonlinear SDE over the interval [0, t], where t 2 ½0; T. Then, operational matrices of integration (omi) based on Bernoulli polynomials are utilized to significantly reduce the complexity of solving obtained SIE through converting such SIE into a nonlinear system of algebraic equations. Error analysis and convergence of suggested technique are also analyzed under some mild conditions. It is concluded that by increasing n^, the approximate solution more accurately estimates the exact solution, where n^denotes the number of elements in the used Bernoulli vector. Finally, some test problems are included to emphasize that the introduced idea is accurate, efficient and applicable. One of the most important innovations of this paper is the numerical simulation of vofBm, which is done in two steps. In the first step, standard Brownian motion (sBm) is constructed via spline interpolation scheme and then block-pulse and hat functions are used in the second step to simulate vofBm. Moreover, the obtained numerical results are also compared with achieved results from Chebyshev cardinal wavelets (Ccw) method to confirm the superiority of the presented method respect to the previous methods. This paper finishes with presenting a real application of this model and solving it via our method.
