In this study, a new hybrid method is developed to solve linear or nonlinear systems of fractional differential equations using Bernoulli wavelets (Bws) and the least squares support vector regression (LS‐SVR). The numerical methods based on operational matrices for solving various kinds of fractional equations have been widely studied in the last decade. In contrast to the existing methods, here we derive the operator of fractional integration, aiming to remove the approximation error. For this purpose, we present Bw operator of Riemann–Liouville fractional integration and use it in our scheme. In the proposed technique, we approximate the unknown functions via Bws, and then with the help of Bw operator of fractional integration and LS‐SVR, we reduce the problem to an algebraic system. In this way, we simplify the computation of the considered system. The error analysis of our method is proposed. Finally, we demonstrate the applicability of the present scheme by solving several numerical examples.