We apply the Monte Carlo method to solving the Dirichlet problem of linear parabolic equations with fractional
Laplacian. This method exploits the idea of weak approximation of related stochastic differential equations
driven by the symmetric stable Lévy process with jumps. We utilize the jump-adapted scheme to approximate
Lévy process which gives exact exit time to the boundary. When the solution has low regularity, we establish
a numerical scheme by removing the small jumps of the Lévy process and then show the convergence order. When
the solution has higher regularity, we build up a higher-order numerical scheme by replacing small jumps with a
simple process and then display the higher convergence order. Finally, numerical experiments including
ten- and one hundred-dimensional cases are presented, which confirm the theoretical estimates and show the numerical
efficiency of the proposed schemes for high-dimensional parabolic equations.