Patch antennas on ferrite substrates allow for pattern control, frequency shifting, and scattering reduction. This is achieved by external magnetic field biasing coupled with the inherent magnetization of the ferrite substrate. Measurements and analytical studies based on the method of moments (MoM) have verified these attractive properties of ferrite substrates. However, verification of the analysis is difficult and, furthermore, previous models have relied on uniform biasing across the substrate. We present a hybrid finite element-boundary integral (FE-BI) method, which permits modeling of the true nonuniform bias fields within the substrate for a more accurate prediction of the ferrite patch performance. After validation of the proposed simulation and a demonstration of the inherent properties of the ferrite patch, it is shown that nonuniform biasing is responsible for additional frequency shifts. We also identify the poor condition of the resulting matrix systems and relate this situation to the predictable occurrence of nonpropagating substrate modes. A more robust iterative solver with preconditioning is, therefore, proposed and applied to handle these situations.