Magnetic surveying encounters challenges in high-precision detection and inversion on strongly magnetized bodies. Non-uniform magnetization caused by self-demagnetization and mutual magnetization makes magnetic field calculation using analytical formulas impractical. The previous discrete numerical methods have suffered from computational inefficiency or low accuracy. To resolve this issue, we develop an algorithm to solve the magnetic field integral equation for calculating the magnetization state of high-susceptibility bodies, combining the improved spatial convolution forward approach with the adaptive relaxation iteration method. To address the excessive computational burden in the magnetization field between 3-D voxel arrays of the discrete model, we construct the circular convolution kernel matrix directly when the subsurface space is divided with an equidistant grid, significantly reducing redundant computations. Then, the fast Fourier transform is used to accelerate the forward discrete circular convolution operation. In addition, we derive an iterative computation scheme that adaptively adjusts the relaxation factor based on magnetic field changes to correct the effective magnetization for algorithm convergence. From a pragmatic perspective, equating the cuboid to equal-volume sphere subdivisions is proposed to enhance computational efficiency further by about 50 per cent despite sacrificing slight accuracy. The sphere and spherical shell models validate the algorithm accuracy, efficiency, and convergence. Furthermore, we analyze the characteristics of the mutual magnetization effect among magnetic bodies under different magnetization conditions. Finally, we demonstrate the applicability of the algorithm through a realistic example. The proposed algorithm shows promise for precisely exploring highly magnetized minerals, underground pipelines, and unexploded ordnance.