In this paper we build the fractional hybrid function of block-pulse functions and Fibonacci polynomials (FHBPF) to numerically solve a class of multiterm variable-order fractional differential equations. We consider fractional derivatives in the Caputo sense and fractional integrals in the Riemann-Liouville sense. We construct an exact integral operator attached to the FHBPF, based on the incomplete beta functions. With the help of the Newton-Cotes collocation procedure, we transform differential equations into systems of algebraic equations, which can be solved by traditional methods such as Newton’s iterative method. We also propose an error investigation method. We conclude the paper with some numerical examples to demonstrate the applicability of the method, but also its simplicity and accuracy.