Using integral equation method, the numerical algorithm of 3D magnetotelluric responses in homogeneous half-space is realized. While calculating tensor Green~ functions , the algorithm of secondary division is adopted to solve the singular problem. The Gaussian quadrature and continued fraction is adopted to calculate the integral with Bessel function. This method not only can ensure the right calculation of tensor Green＇s functions but also can get high precision. The correctness of algorithm has been proved by model comparing and analyse. This efficient algorithm will help to study the distributing laws of 3D magenetotelluric responses and to setup basis for reaserch of 3D inversion.