The tomographic AP method is so far the best method in separating the Alcock-Paczynski (AP) signal from the redshift space distortion (RSD) effects and deriving powerful constraints on cosmological parameters using the $lesssim40h^{-1} rm Mpc$ clustering region. To guarantee that the method can be easily applied to the future large scale structure (LSS) surveys, we study the possibility of estimating the systematics of the method using fast simulation method. The major contribution of the systematics comes from the non-zero redshift evolution of the RSD effects, which is quantified by $hatxi_{Delta s}(mu,z)$ in our analysis, and estimated using the BigMultidark exact N-body simulation and approximate COLA simulation samples. We find about 5%/10% evolution when comparing the $hatxi_{Delta s}(mu,z)$ measured as $z=0.5$/$z=1$ to the measurements at $z=0$. We checked the inaccuracy in the 2pCFs computed using COLA, and find it 5-10 times smaller than the intrinsic systematics of the tomographic AP method, indicating that using COLA to estimate the systematics is good enough. Finally, we test the effect of halo bias, and find $lesssim$1.5% change in $hatxi_{Delta s}$ when varying the halo mass within the range of $2times 10^{12}$ to $10^{14}$ $M_{odot}$. We will perform more studies to achieve an accurate and efficient estimation of the systematics in redshift range of $z=0-1.5$.