The paper is devoted to the numerical solutions of fractional PDEs based on its probabilistic interpretation, that is, we construct approximate solutions via certain Monte Carlo simulations. The main results represent the upper bound of errors between the exact solution and the Monte Carlo approximation, the estimate of the fluctuation via the appropriate central limit theorem (CLT) and the construction of confidence intervals. Moreover, we provide rates of convergence in the CLT via Berry-Esseen type bounds. Concrete numerical computations and illustrations are included.