Prop. \[prop:lip\_U\] and with the same notation, we have $U^n= \Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k U^k\right)$, ${\tilde U}^n= \Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k {\tilde U}^k+a^n\right)$. Hence $$\begin{aligned}
&\|{\tilde U}^n-U^n\|_\infty= \left\|\Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k {\tilde U}^k+\rho_\alpha a^n\right)
-\Psi\left(\sum_{k=n+1}^{N} {\overline c}_{n}^k  U^k\right)\right\|_\infty\\
&\le \sum_{k=n+1}^{N}{\overline c}_{n}^k \left\| {\tilde U}^k - U^k\right\|_\infty+{\rho_\alpha}\|a^n\|_\infty.\end{aligned}$$ The previous inequality is equivalent to $$\bar D_{\Delta t}^{\alpha}(\|{\tilde U}^n-U^n\|_\infty)\le \|a^n\|_\infty.$$ We conclude, by the backward version of discrete Gronwall inequality , the uniform convergence of $u_{\sigma}$ to $u$. The convergence of the gradients follows by .

In this paper we have focused on the numerical approximation of the classical solution of the PDE system and we have assumed the existence of a unique classical solution. We refer to [@q] for the existence and uniqueness of a weak solution to the system, where the solution of the fractional Fokker-Planck equation is defined in a time-fractional Bochner space and the equation is considered in the sense of distributions. It would be very interesting to define some notion of renormalized solution to fractional PDEs and consider weak solution to the system, as well as approximation schemes. In the standard case, convergence of a finite differences scheme to the weak solution of a Mean Field Games system was shown in 