  A random walk-based method is proposed to efficiently compute the solution of
a large class of fractional in time linear systems of differential equations
(linear F-ODE systems), along with the derivatives with respect to the system
parameters. Such a method is unbiased and unconditionally stable, and can
therefore be used to provide an unbiased estimation of individual entries of
the solution, or the full solution. By using stochastic differentiation
techniques, it can be used as well to provide unbiased estimators of the
sensitivities of the solution with respect to the problem parameters without
any additional computational cost. The time complexity of the algorithm is
discussed here, along with suitable variance bounds, which prove in practice
the convergence of the algorithm. Finally, several test cases were run to
assess the validity of the algorithm.
