This paper presents a mathematical model for the 3D nonlinear magnetostatic field based on integral equations with fictitious surface and volume magnetic charges. The solution is performed by the Extended Boundary Element Method including Surface Elements (SE) and Volume Elements (VE). Examples of the calculation for both linear and nonlinear magnetic systems are presented.