In this paper, we formulate a numerical method to approximate the solution of two‐dimensional optimal control problem with a fractional parabolic partial differential equation (PDE) constraint in the Caputo type. First, the optimal conditions of the optimal control problems are derived. Then, we discretize the spatial derivatives and time derivatives terms in the optimal conditions by using shifted discrete Legendre polynomials and collocations method. The main idea is simplifying the optimal conditions to a system of algebraic equations. In fact, the main privilege of this new type of discretization is that the numerical solution is directly and globally obtained by solving one efficient algebraic system rather than step‐by‐step process which avoids accumulation and propagation of error. Several examples are tested and numerical results show a good agreement between exact and approximate solutions.