## 2 A New Class of Linear Multistep Formulae  
### Construction of a New Class of Linear Multistep Formulae by Classic Differential Quadrature  
The initial-value problem of an ordinary differential equation was shown in the following:  
\[\left\{\begin{array}{l}\dot{x}(t)=f(t,x)\\ x(t=0)=x_{0}\end{array}\right. \tag{1}\]  
where \(t\) is time and \(x\) is a one-dimensional variable; \(f(t,x)\) is a function of the time \(t\) and the variable \(x\); and \(x_{0}\) is the initial value of an unknown state variable.  
One of main methods for solving the initial-value problem, as indicated in Equation (1), is the use of classical linear multistep formulae (LMF) [11], and one kind of the LMF is the following backward differentiation formula (BDF), which was described as following:  
\[\sum_{i=0}^{s}\alpha_{i}x_{n+i}=hf_{n+s} \tag{2}\]  
where \(s\) is the number of steps; \(\alpha_{i}\) is integration coefficient; \(h\) denotes the integration step; \(x_{n+i}\approx x(t_{n+i})\) is the approximate value of a state variable at a time point \(t_{n+i}\); and \(f_{n+s}=f(t_{n+s},x_{n+s})\).  
Obviously, when \(s\) takes one, Equation (2) is the first-order backward Euler method that is A-stable and L-stable. Based on the traditional BDF, new linear multistep methods were constructed as follows:  
\[\sum_{i=0}^{s}\widetilde{\alpha}_{i}x_{n+i}=h\beta_{s}f_{n+s} \tag{3}\]  
where \(\widetilde{\alpha}_{i}\) and \(\beta_{s}\) are integration coefficients of new LMF, respectively.  
In order to obtain the unknown integration coefficients of new LMF, we similarly used \(s\)-stage \(s\)-order time-domain differential quadrature method (DQM) [12; 13; 14] to solve Equation (1). For any time-step (\(h=t_{n+1}-t_{n}\)), we have the following equation:  
\[G\left[\begin{array}{c}\widetilde{x}_{n+1}\\ \vdots\\ \widetilde{x}_{n+s}\end{array}\right]+G_{0}x_{n}=h\left[\begin{array}{c}f \big{(}\widetilde{t}_{1},\widetilde{x}_{n+1}\big{)}\\ \vdots\\ f\big{(}\widetilde{t}_{s},\widetilde{x}_{n+s}\big{)}\end{array}\right] \tag{4}\]  
where \(c_{i}=(\widetilde{t}_{i}-t_{n})/h,i\in(0,s)\) is the non-equidistant grid point, which is distributed over the regularization interval \([0,1]\); \(\widetilde{t}_{0}=t_{n}\) is a time-starting point; \(\widetilde{t}_{s}=t_{n+1}\) is a time-end point; \(\widetilde{t}_{i}=t_{n}+c_{i}h,i\in(1,s)\) are internal points; and \(\widetilde{x}_{n+i}\approx x(\widetilde{t}_{i})\) is the approximate value of a state variable at a time-internal point \(\widetilde{t}_{i}\).  
\[G=\left[\begin{array}{ccc}\delta_{11}&\cdots&\delta_{1s}\\ \vdots&\ddots&\vdots\\ g_{s1}&\cdots&g_{ss}\end{array}\right],G_{0}=\left[\begin{array}{c}\delta_{1 0}\\ \vdots\\ g_{s0}\end{array}\right] \tag{5}\]where \(G\) is called the weighted coefficient matrix of DQM, and the calculation formulae of each element in \(G\) is given in Reference [15]\(G_{0}=-Ge\), where \(e\) is the column vector with each element equaling to one.  
In Equation (4), the grid point distribution has to be a non-equidistant grid, including Chebyshev grid, Chebyshev-Gauss-Lobatto grid, and Legendre grid [16]. In this paper, we used the DQM with Chebyshev grid points [16] for constructing the new LMF. Suppose:  
\[\widetilde{\alpha}_{i}=g_{si}/g_{ss},i\in(0,s),\beta_{s}=\sum_{i=1}^{s}i \widetilde{\alpha}_{i} \tag{6}\]  
where \(\widetilde{\alpha}_{s}=1\).  
As an example, the new three-step LMF using Chebyshev grid points was shown below:  
\[\frac{15-4\sqrt{2}}{9}hf_{n+3}=-\frac{1}{9}x_{n}+\frac{4\sqrt{2}-4}{9}x_{n+1} -\frac{4\sqrt{2}+4}{9}x_{n+2}+x_{n+3} \tag{7}\]  
The new four-step LMF using Chebyshev grid points is as follows:  
\[\frac{108-32\sqrt{3}}{57}hf_{n+4}=\frac{1}{19}x_{n}+\frac{16\sqrt{3}-32}{57}x _{n+1}+\frac{4}{57}x_{n+2}-\frac{16\sqrt{3}+32}{57}x_{n+3}+x_{n+4} \tag{8}\]  
In addition, the backward Euler method is obtained when \(s\) takes one in Equation (3) regardless of the kind of the used grid. Similarly, when \(s\) takes two, we get the second-order Gear method.