This thesis describes a method to integrate the SPICE device models in Obreshkov based circuit simulation. Recently, a high-order A-stable and L-stable Obreshkov based method was proposed to simulate the transient response for general nonlinear circuits, which provides more than one order of magnitude speedup as compared to the traditional second order numerical integration methods used in the SPICE engine, without sacrificing the accuracy. This method requires calculation of higher order derivatives of Nodal equations, with orders higher than 2. This thesis describes the process to build a rooted tree from the source code of SPICE device models, and how to use this rooted tree to calculate high-order derivatives needed in the Obreshkov formula. This method was tested on the SPICE diode model, and the calculated values of the derivatives have been found to match the expected values up to 6 decimal places.