This study investigates differential equations characterized by multiple fractional derivatives defined in the Caputo sense. To address this problem, we employ a finite difference method, approximating the fractional derivatives using the L1 and L2 schemes on a uniform mesh. A rigorous error analysis is conducted to derive theoretical estimates of the numerical scheme&rsquos accuracy. Numerical computations are performed to evaluate the scheme&rsquos performance and precision, revealing that the rate of convergence achieves a linear order, consistent with the established theoretical estimates. To enhance the accuracy of the scheme, we implement Richardson extrapolation, demonstrating an effective strategy for improving the method&rsquos accuracy.