The main purpose of this work is to develop effective spectral collocation methods for solving numerically general form of non-linear fractional two-point boundary value problems with two end-point singularities. Specifically, we define two classes of Müntz-Legendre polynomials to deal with the singular behaviors of the solutions at both end-points, which enhance greatly the accuracy of the numerical solutions. Important and practical formulas for ordinary derivatives of Müntz-Legendre polynomials are derived. Then using a three-term recurrence formula and Jacobi-Gauss quadrature rules, applicable and stable numerical approaches for computing left and right Caputo fractional derivatives of these basis functions are provided. One of the important features of the present method is to show how to recover spectral accuracy from the end-point singularities for the coupled systems of fractional differential equations with left and right fractional derivatives using the mentioned approximation basis. Moreover, we present a detailed analysis with accurate convergence rates of the singular behavior of the solution to a rather general system of nonlinear fractional differential equations including left and right fractional derivatives. Numerical results show the expected convergence rates. Finally, some numerical examples are given to demonstrate the performance of the proposed schemes.