We introduce a numerical algorithm based on the hybrid of block-pulse functions and fractional-order Jacobi polynomials in order to solve fractional differential equations. The fractional derivative is described in the Caputo sense. The Riemann–Liouville fractional integral operator for these basis functions is constructed. This result together with the shifted Gauss–Chebyshev collocation points are utilized to reduce the original problem to a system of nonlinear algebraic equations. By means of solving the given system, the numerical solution of the main problem is derived. Then, the method is applied for solving the Bagley–Torvik initial and boundary value problems. After that, an error estimation is presented for the expansion of a given function based on the fractional-order basis functions. Finally, for demonstrating the precision and good performance of the new method, several numerical examples are considered and the results are compared with the exact or approximate solutions obtained by other existing techniques.
