Accurate magnetic field calculation is the premise of electromagnetic performance prediction. Conventional subdomain (SD) techniques assume that the iron’s relative permeability is infinite, leading to falsely overestimated flux density. We propose an accurate magnetic field analytical model for permanent magnet (PM) in-wheel machines considering iron’s magnetization nonlinearity and saturation. Specifically, according to the excitation source and topology, the entire solution domain of the machine is divided into sub-regions such as stator slots/teeth, stator slot-openings/tooth-tips, air-gap, and rotor slots/teeth, etc. Poisson’s or Laplace’s magnetic vector potential (MVP) equations are solved using Maxwell’s electromagnetic theory and complex Fourier series methods in each sub-region. Specifically, in our approach, The Cauchy product theorem addresses the discontinuous magnetic permeability change at the slot and tooth interface. The machine’s magnetic saturation effect is considered by combining the actual magnetization characteristics of iron with an iterative algorithm. The general solution for the MVP is solved using the boundary conditions between adjacent subregions. Subsequently, electromagnetic properties such as air-gap flux density, back electromotive force (EMF), and electromagnetic torque are obtained. The accuracy of the analytical model is verified by finite element analysis (FEA) and prototype tests, which proved that the proposed analytical model can consider the iron’s nonlinearity and the magnetic saturation. In addition, the inaccurate overestimation of electromagnetic torque and air-gap magnetic flux density by the conventional SD techniques has also been proven.