  The equation with the time fractional substantial derivative and space
fractional derivative describes the distribution of the functionals of the
L\'evy flights; and the equation is derived as the macroscopic limit of the
continuous time random walk in unbounded domain and the L\'evy flights have
divergent second order moments. However, in more practical problems, the
physical domain is bounded and the involved observables have finite moments.
Then the modified equation can be derived by tempering the probability of large
jump length of the L\'evy flights and the corresponding tempered space
fractional derivative is introduced. This paper focuses on providing the high
order algorithms for the modified equation, i.e., the equation with the time
fractional substantial derivative and space tempered fractional derivative.
More concretely, the contributions of this paper are as follows: 1. the
detailed numerical stability analysis and error estimates of the schemes with
first order accuracy in time and second order in space are given in
{\textsl{complex}} space, which is necessary since the inverse Fourier
transform needs to be made for getting the distribution of the functionals
after solving the equation; 2. we further propose the schemes with high order
accuracy in both time and space, and the techniques of treating the issue of
keeping the high order accuracy of the schemes for {\textsl{nonhomogeneous}}
boundary/initial conditions are introduced; 3. the multigrid methods are
effectively used to solve the obtained algebraic equations which still have the
Toeplitz structure; 4. we perform extensive numerical experiments, including
verifying the high convergence orders, simulating the physical system which
needs to numerically make the inverse Fourier transform to the numerical
solutions of the equation.
