A new numerical integral method of the Hankel transform which included the dual Bessel function in the integrand was presented. Because of the strong-oscillation and slow-decay properties of the Bessel functions, it is impossible to use a normal numerical integration method. In our method, the integral interval [0, ∞) is divided into two parts: [0, λ
_0
] and [λ
_0
, ∞). In [0,λ
_0
], the integration can be calculated using a Gauss-Legendre quadrature method between the adjacent zero points. And in [λ
_0
, ∞), the Hankel transform is converted to the Fourier cosine or sine transform, and can be calculated with a polygonal approximation algorithm. The proper value of the dividing point ofλ
_0
was discussed. The simulation results of the test functions show that the algorithm is correct and accurate. Finally, the forward problem of the multi-layer medium magnetic induction tomography was calculated with this new method, and the eddy current distribution in the model and the induced voltages in the detection coils were simulated.
