In this article, we consider the numerical solution for Poisson's equation in axisymmetric geometry. When the boundary condition and source term are axisymmetric, the problem reduces to solving Poisson's equation in cylindrical coordinates in the two-dimensional ($r,z$) region of the original three-dimensional domain $S$. Hence, the original boundary value problem is reduced to a two-dimensional one. To make use of the Mechanical quadrature method (MQM), it is necessary to calculate a particular solution, which can be subtracted off, so that MQM can be used to solve the resulting Laplace problem, which possesses high accuracy order $\mathcal{O}$($h^3_{max}$) and low computing complexities. Moreover, the multivariate asymptotic error expansion of MQM accompanied with $\mathcal{O}$($h^3_i$) for all mesh widths $h_i$ is got. Hence, once discrete equations with coarse meshes are solved in parallel, the higher accuracy order of numerical approximations can be at least $\mathcal{O}$($h^5_{max}$) by the splitting extrapolation algorithm (SEA). Meanwhile, a posteriori asymptotic error estimate is derived, which can be used to construct self-adaptive algorithms. The numerical examples support our theoretical analysis.