  In this paper, we study the numerical approximation of a system of PDEs with
fractional time derivatives. This system is derived from an optimal control
problem for a time-fractional Fokker-Planck equation with time dependent drift
by convex duality argument. The system is composed by a time-fractional
backward Hamilton-Jacobi-Bellman and a forward Fokker-Planck equation and can
be used to describe the evolution of probability density of particles trapped
in anomalous diffusion regimes. We approximate Caputo derivatives in the system
by means of L1 schemes and the Hamiltonian by finite differences. The scheme
for the Fokker-Planck equation is constructed such that the duality structure
of the PDE system is preserved on the discrete level. We prove the well
posedness of the scheme and the convergence to the solution of the continuous
problem.
