The trend of higher integration density in VLSI has made the time-domain circuit simulation a bottleneck in today's IC design flows. For designs with millions of elements, SPICE-like simulation can easily take days, even weeks, to complete the task. Therefore, accurate yet fast circuit simulation methods have always been one of the major demands in industry. A different aspect of solving the simulation problem is proposed in this dissertation. The proposed method, called Matrix Exponential Method (MEXP), is based on the analytical solution of the ordinary different equations. The MEXP method has the benefits of accuracy, stability, scalability, parallelizability and adaptivity due to its analytical nature and simple kernel computation. To utilize the adaptivity of MEXP, an algorithm is proposed to dynamically determine locally optimal step size and order of Krylov subspace at each time step to reduce the overall computation in the simulation. Moreover, the scaling effect of MEXP is observed. The scaling effect enables MEXP to simulate a circuit with larger step size as the time instant moves forward. With the scaling effect, MEXP demonstrates the comparable capability of stiffness handling as the widely adopted backward Euler and trapezoidal methods. The parallelism of MEXP is also demonstrated. The GPU implementation is shown in this dissertation. The restarted MEXP is also proposed for GPU implementation to mitigate the memory usage on resource restricted environment. Finally, the rational Krylov subspace method is investigated. Instead of using A as the basis, the rational Krylov subspace utilizes (I-[gamma]A)⁻¹ as the basis. The rational basis transforms the spectrum of a matrix to relax the stiffness constraint in the Krylov subspace method. Furthermore, the rational Krylov subspace can avoid the regularization process of singular C. This dissertation applies MEXP with the rational Krylov subspace on the power distribution network (PDN) design. Combining with the capability of adaptivity and accuracy of MEXP, MEXP with the rational Krylov subspace could exploit a large step size in the low-frequency response to further improve the performance of the power grid simulation