We present a novel parallel algorithm to solve fractional-order systems involving Caputo-types derivative. The numerical method used for finding the solution is Adams-Bashforth-Moulton (ABM) predictor-corrector scheme. Using MPI and OPENMP frameworks the algorithm was implemented to run on a BlueGene/P supercomputer. The same algorithm that ran on the BlueGene/P cluster was adapted to run on GPU and exploit the GPU's capabilities. A comparison between these approaches is presented by running numerical simulations of the same three-dimensional fractional-order systems.

The algorithm implemented in CUDA that approximates the solution using the ABM method was also adapted to use Diethelm's method. Having the same algorithm that runs on the same system, enables an accurate and practical comparison between the numerical methods.