The aim of this paper is to investigate the numerical solutions of delayed fractional differential equations (FDEs). A constant time delay is taken in the Caputo derivative of the state variable and also in the state variable. To obtain solutions, we extend the block boundary value method (BBVM) and show the convergence of the resulting scheme. Further, it is found that the order of convergence is $\min\{m,q-\delta+1\}$, where $m$ and $q$ are the order and block-size of the BBVM, respectively, and $\delta$ lies in between $0$ and $1$. In our methodology, we approximate the fractional-ordered derivative by a combination of the $m^{\text{th}}$-order BBVM and a $q^{\text{th}}$-order Lagrange interpolating polynomial. Lastly, we analyze the global stability of the numerical scheme and with the aid of some numerical examples, its computational effectiveness is proved. In examples, we analyse some fractional order delay differential systems and present a comparison of solutions for different fractional order derivatives. Furthermore, it is demonstrated that the global errors diminish as the order of the fractional derivative decreases, and this fact is supported by the theoretical order of convergence.
MSC Classification: 34K37 , 34D23 , 65L20