  Many differential equations with physical backgrounds are described as
gradient systems, which are evolution equations driven by the gradient of some
functionals, and such problems have energy conservation or dissipation
properties. For numerical computation of gradient systems, numerical schemes
that inherit the energy structure of the equation play important roles, which
are called structure-preserving. The discrete gradient method is one of the
most classical framework of structure-preserving methods, which is at most
second order accurate. In this paper, we develop a higher-order
structure-preserving numerical method for gradient systems, which includes the
discrete gradient method. We reformulate the gradient system as a coupled
system and then apply the discontinuous Galerkin time-stepping method.
Numerical examples suggests that the order of accuracy of our scheme is $(k+1)$
in general and $(2k+1)$ at nodal times, where $k$ is the degree of polynomials.
