Abstract: A new computational method to solve the hyperbolic (gyrokinetic Vlasov) equation and the elliptic (Poisson-like) equation at the polar axis is proposed. It is shown that the value of a scalar function at the polar axis can be predicted by its neighbouring values based on the continuity condition. This continuity condition systematically solves the pole problems including the singular factor 1/r in the hyperbolic equation and the inner boundary in the elliptic equation. The proposed method is applied to the global gyrokinetic simulation of the tokamak plasma with the magnetic axis included.