To analyze the vibration characteristics of transformer windings with high accuracy, a solution approach combining the axial dynamic model of windings with the two-stage fourth-order Gauss implicit Runge Kutta （GRK24） method is developed. The simplified model of "mass-spring-damping" is adopted for the axial dynamics modeling of windings, and the vibration equation is solved by the matrix form GRK24 （GRK24-M） method in which the displacement, velocity and acceleration solutions can reach fourth-order accuracy simultaneously. The windings structure of a 50 kVA transformer is taken as an example. The effectiveness of the proposed approach solving the windings vibration responses under the impact load and the electromagnetic force load is verified. The superiority of the proposed GRK24-M method is also illustrated by comparing the solution results with the Newmark method and the HHT-α method.