Numerical instabilities in the magnetic field closed-form representation of homogeneous cuboid permanent magnets

We investigate the behavior of the well-known closed-form representation of the magnetic field of a cuboid permanent magnet with homogeneous magnetization derived in [1, 2]. These formulas are of practical value, e.g. for the design of magnet systems, experiment layout [2], state computation in sensor systems [3], and they are also implemented in the magnet and sensor system design open-source package Magpylib [4].
However, there is little awareness for the numerical limitations, when evaluating these formulas in 64bit floating-point representation along the cuboid edges and their spatial extensions to infinity. While the analytical representation becomes an indeterminate form along those lines, which is easily computed in the limit, a floating-point evaluation already breaks down in the vicinity of these special cases.
This is demonstrated in the figures below, where some of the instabilities have been removed via symmetry transformation to estimate the error. The size of the numerical error is threatening when applications like NMR require a relative precision of 10-6.

[1] J. M. Camacho and V. Sosa, Rev. Mex. Fıs. E 59, 8 (2013).
[2] Z. J. Yang, T. Johansen, H. Bratsberg, G. Helgesen and A. T. Skjeltorp, Supercond. Sci. Technol. 3, 591 (1990).
[3] D. Cichon, R. Psiuk, H. Brauer and H. Topfer, IEEE Sens. J. 19, 2509 (2019).
[4] M. Ortner and L. G. Coliado Bandeira, SoftwareX 11, 100466 (2020); https://magpylib.readthedocs.io/en/latest/