A computational methodology, based on the coupling of the finite element and boundary element methods, is developed for the solution of magnetothermal problems. The finite element formulation and boundary element formulation, along with their coupling, are discussed. The coupling procedure is also presented, which entails the application of the LU decomposition to eliminate the need for the direct inversion of matrices resulting from FE‐BE formulation, thereby saving computation time and storage space. Corners for both FE‐BE interface and BE regions, where discontinuous fluxes exist, are treated using the double flux concept. Numerical results are presented for three different systems and compared with analytical solutions when available. Numerical experiments suggest that for magnetothermal problems involving small skin depths, a careful mesh distribution is critical for accurate prediction of the field variables of interest. It is found that the accuracy of the temperature distribution is strongly dependent upon that of the magnetic vector potential. A small error in the magnetic vector potential can produce significant errors in the subsequent temperature calculations. Thus, particular attention must be paid to the design of a suitable mesh for the accurate prediction of vector potentials. From all the cases examined, 4‐node linear elements with adequate progressive coarsening of meshes from the surface gave the results with best accuracy.