Discontinuities in time-based simulation

I am trying to simulate a fault current limiter (saturable reactor) in a test circuit. When a fault occurs, the reactor voltage has discontinuities that correspond to locations in time when the current transitions to/from nominal levels to faults levels (reactor moves in and out of saturation).

Previous authors show these discontinuties as small, but continous humps at the voltage waveform peaks. I'm becoming more convinced that this is a numerical issue. Could it be possible that a more ruggidized solver would help to reduce the influence of the doscontinuities? The B/H curve and the dBdH vs H curves are smooth.

I've used both sequential and trapezoidal integration technigues to perform the simulation. Any suggestions? I have plots if there is a way to share them.

Thanks :) M