This study introduces a novel numerical method for approximating variable-order fractional Brownian motion, representing a significant advancement in the field of stochastic processes. The proposed method enhances the modeling accuracy of complex phenomena by accommodating variable-order Brownian motion. Additionally, it mitigates the computational challenges typically associated with modeling such processes. The innovative approach employs a newly developed and straightforward matrix-based algorithm grounded in B-spline functions, offering an efficient, accurate, and computationally simple technique for approximating variable- order fractional Brownian motion. Also, this study focuses on solving a novel class of integral equations driven by variable-order fractional Brownian motion. The proposed method uses the features of barycentric rational interpolants and the spectral method to provide a simple and accurate approach, thereby reducing the complexities associated with solving such integral equations. The convergence of the method is analyzed in detail, and its theoretical robustness is emphasized. Furthermore, several numerical experiments have been conducted, demonstrating the reliability and adaptability of the method in challenging stochastic models. All numerical results have been analyzed using statistical methods to ensure greater reliability and accuracy.