In your MOSFET model, the parameter “Ld” is much too big. If it’s greater then half the channel length “L”, the effective length becomes negative, which is, of course, unphysical.
Anyway, the problem here is that the trapezoidal method can’t handle steps quite well. Please choose the integration method GEAR3. I hope that works better.