<p>This work introduces a numerical method based on modified cubic B-splines for numerically solving a fractional sub-diffusion model. We analyse the considered model with one parameter Mittag-Leffler function as a kernel that is non-local in nature and non-singular. The proposed method considers a modified cubic B-spline collocation approach in space and an efficient finite difference scheme in time. The unconditional stability of the method is established through the von Neumann technique. Further, the method is proved to be second order convergent in both space and time. Numerical results support the theory and show the applicability of the proposed method.</p>