  We introduce a new numerical method, based on Bernoulli polynomials, for
solving multiterm variable-order fractional differential equations. The
variable-order fractional derivative was considered in the Caputo sense, while
the Riemann--Liouville integral operator was used to give approximations for
the unknown function and its variable-order derivatives. An operational matrix
of variable-order fractional integration was introduced for the Bernoulli
functions. By assuming that the solution of the problem is sufficiently smooth,
we approximated a given order of its derivative using Bernoulli polynomials.
Then, we used the introduced operational matrix to find some approximations for
the unknown function and its derivatives. Using these approximations and some
collocation points, the problem was reduced to the solution of a system of
nonlinear algebraic equations. An error estimate is given for the approximate
solution obtained by the proposed method. Finally, five illustrative examples
were considered to demonstrate the applicability and high accuracy of the
proposed technique, comparing our results with the ones obtained by existing
methods in the literature and making clear the novelty of the work. The
numerical results showed that the new method is efficient, giving high-accuracy
approximate solutions even with a small number of basis functions and when the
solution to the problem is not infinitely differentiable, providing better
results and a smaller number of basis functions when compared to
state-of-the-art methods.
