This paper presents the development of a fractional hybrid function composed of block-pulse functions and Fibonacci polynomials (FHBPF) for the numerical solution of multiterm variable-order fractional differential equations. By replacing \(x\rightarrow x^{\alpha }\) in FHBPF and utilizing incomplete beta functions, we construct the method with a focus on fractional derivatives in the Caputo sense and fractional integrals in the Riemann–Liouville sense. A key advantage of the proposed method is its exact computation of the Riemann–Liouville fractional integral operator. Using the Newton–Cotes collocation method, we transform the differential equations into systems of algebraic equations, which are then solved with traditional techniques such as Newton’s iterative method. An error analysis method is also introduced. Several numerical examples are provided to demonstrate the effectiveness and simplicity of the approach, highlighting its efficiency even for relatively small base sizes.
