In this paper, we introduce a method based on Radial Basis Functions (RBFs) for the numerical approximation of Mathieu differential equation with two fractional derivatives in the Caputo sense. For this, we suggest a numerical integration method for approximating the improper integrals with a singularity point at the right end of the integration domain, which appear in the fractional computations. We study numerically the affects of characteristic parameters and damping factor on the behavior of solution for fractional Mathieu differential equation. Some examples are presented to illustrate applicability and accuracy of the proposed method. The fractional derivatives order and the parameters of the Mathieu equation are changed to study the convergency of the numerical solutions.