In this paper, we propose an implicit-explicit Monte Carlo method (IMEX-MCM) for solving semi-linear parabolic partial differential equations (PDEs) driven by
α
\alpha
-stable Lévy process on bounded domains, also known as semi-linear PDEs involving integral fractional Laplacian operator (i.e.,
α
∈
(
0
,
2
)
\alpha \in (0,2)
) and classical Laplacian operator (i.e.,
α
=
2
\alpha =2
). The main idea behind our approach is to divide the time interval into
N
N
sub-intervals and derive the Feynman-Kac formula within each time sub-interval. Then, we develop a fully discretized parallel scheme using the implicit-explicit technique for temporal discretization, followed by applying a novel walk-on-sphere method to obtain the numerical solution. Notably, the proposed method takes advantage of the Poisson kernel associated with the stochastic processes to accommodate the full discretization of the classical and fractional semi-linear PDEs in a single framework. Moreover, we rigorously analyze the error bound for the proposed scheme that is totally explicit with respect to the number of simulation paths
M
M
and the time step size
Δ
t
\Delta t
. As a by-product, we prove the maximum principle of the proposed schemes for fractional Allen Cahn equation with double well potential and provide an error estimate by using the established maximum principle. Extensive numerical experiments in two and three-dimensional spaces are presented to verify our theoretical results and to show the robustness and accuracy of our method.