In this study, we present a practical numerical approach for solving a fractional-order diffusion problem, and extend it to address the optimality system of a fractional-order diffusion problem. The methodology involves numerical analysis of the solution to the optimal control problem within a fractional diffusion system. We utilize a finite difference method on a bounded domain, considering the fractional time derivative in a Riemann–Liouville sense. The discretisation of both state and adjoint equations forms the basis for developing numerical algorithms. The obtained results are then analysed, including an examination of convergence properties and stability under different conditions. To illustrate the applicability of our approach, we provide a numerical example. This example serves as a practical demonstration, showcasing the capabilities and insights offered by our numerical scheme.