In this paper a modified model of ferromagnetic hysteresis suitable for mixed-signal simulations in VHDLAMS is presented. The aim of this paper is to demonstrate how a numerically stable and accurate implementation of the Jiles-Atherton model can be achieved using a 4th order Runga-Kutta integration of the derivative of magnetization with respect to the field strength (H). While most SPICE-like implementations require inconvenient integration in time to obtain the magnetization derivative, our approach is more general as it does not rely on the underlying differential equation solver for this purpose. The model addresses the non-physical situation of negative BH slopes and proposes an alternative implementation of the anhysteretic function using a polynomial approximation of the Langevin function for low signal levels and a new function with no discontinuities. Model efficiency is improved by monitoring the change in H and only activating the integration function when H changes by a specified amount.