<p>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) with an inhomogeneous fractional exponent, 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 also provide unbiased estimators of the sensitivities of the solution with respect to the problem parameters. This is achieved with only a minimal increase in computational cost and negligible memory overhead, since the sensitivities are computed within the same random walk framework without requiring additional independent simulations or separate solvers. 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 with random matrices were run to assess the validity of the algorithm, together with a more realistic example.</p>