  A new computational method to solve the hyperbolic (Vlasov) equation coupled
to 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.
