We consider the discrete gradient method for dissipative linear-gradient systems, which strictly replicates the dissipation property, yielding a remarkable stability. However, it also replicates the nonlinearity of an original equation. To overcome this, we can employ multistep linearly implicit schemes as a relaxation; however, it can in turn destroy the originally aimed stability. Matsuo-Furihata (2014) introduced a dynamical systems viewpoint to understand the behavior for a toy scalar problem. In this letter, we show that their method can work also for the two-dimensional Duffing equation. There a new concept of semi-strong Lyapunov functionals is required.