•The method solves multi-order and variable coefficients RLFDE and CFD.•Theorem 3.1. transforms RLFI to a linear algebraic equation.•Corollary 3.1.1. transforms CFD to a linear algebraic equation.•Theorem 3.2. allows to find with for RLFDE.•Corollary 3.2.1. allows to find and for CFDE. In this paper, a numerical method is developed to obtain a solution of Caputo’s and Riemann-Liouville’s Fractional Differential Equations (CFDE and RLFDE). Scientific literature review shows that some numerical methods solve CFDE and there is only one paper that numerically solves RLFDE. Nevertheless, their solution is limited or the Fractional Differential Equation (FDE) to be solved is not in the most general form. To be best of the author’s knowledge, the proposed method is presented as the first method that numerically solves RLFDE which includes multi-order fractional derivatives and variable coefficients. The method converts the RLFDE or CFDE to be solved into an algebraic equation. Each Riemann-Liouville’s or Caputo’s Fractional Derivative (RLFD and CFD), derived from the RLFDE or CFDE respectively, is conveniently written as a set of substitution functions and an integral equation. The algebraic equation, the sets of substitution functions and the integral equations are discretized; and then solved using arrays. Some examples are provided for comparing the obtained numerical results with the results of other papers (when available) and exact solutions. It is demonstrated that the method is accurate and easy to implement, being presented as a powerful tool to solve not only FDE but also a wide range of differential and integral equations.