In this paper, we consider numerical solutions of the fractional diffusion equation with the $α$ order time fractional derivative defined in the Caputo-Hadamard sense. A high order time-stepping scheme is constructed, analyzed, and numerically validated. The contribution of the paper is twofold: 1) regularity of the solution to the underlying equation is investigated, 2) a rigorous stability and convergence analysis for the proposed scheme is performed, which shows that the proposed scheme is $3 + α$ order accurate. Several numerical examples are provided to verify the theoretical statement.