The high-order A- and L-stable Obreshkov-based method was recently proposed for simulating the transient response of general nonlinear circuits in the time-domain. This method has consistently resulted in more than one order-of-magnitude speedup compared with the traditional methods used in the SPICE engine. Nonetheless, the formulation of this Obreshkov-based approach assumed that the memory elements in the circuits are characterized by a linear constituent equation. The goal of this paper is to generalize this approach to handle memory elements that are described by nonlinear constituent equations.