The method based on block pulse functions (BPFs) has been proposed to solve different kinds of fractional differential equations (FDEs). However, high accuracy requires considerable BPFs because they are piecewise constant and not so smooth. As a result, it increases the dimension of operational matrix and computational burden. To overcome this deficiency, a novel numerical method is developed to solve fractional differential equations. The method is based upon hybrid of BPFs and Bernstein polynomials (HBBPs), which are piecewise smooth. The HBBPs operational matrix of fractional‐order integral is derived to reduce the FDEs to a system of algebraic equations. Then the numerical solution of the FDEs is obtained through solving the system of algebraic equations. The convergence analysis is conducted for the suggested scheme, and the upper bound of error of the solution is given. Finally, illustrative examples are presented to demonstrate the validity, applicability, and efficiency of the proposed technique in contrast with other approaches.