Electromagnetic transient (EMT) simulation is one of the routine tasks being executed in modern power systems. Typically, the most common EMT simulation method is the trapezoidal integration method, however, sustained numerical oscillations are possible to appear using this method, for instance, getting the voltage waveform across an inductance during a sudden current interruption. Two-stage diagonally implicit Runge-Kutta (2S-DIRK) method is a L-stable numerical method, which is able to avoid numerical oscillation. The 2S-DIRK method has been applied to solve simple RLC circuit simulations, in this paper, it is extended to system level EMT simulations considering the electromagnetic transient characteristics and nonlinearity for generators, transformers and loads. Case studies are performed to verify the effectiveness of the 2S-DIRK method in power system level EMT simulation, by comparing the accuracy and numerical stability with traditional methods.