Page 152 at https://link.springer.com/content/pdf/10.1007%2F978-0-387-73468-2.pdf
Hi readers, I have tried substituting y'(0) , y''(0) , y'''(0) and y''''(0) into equation A that is the first equation. However, still couldn't really get what Maclaurin derived for y(1). Anyone can provide some workings on how Maclaurin did that?
Assuming convergence (so the formula works at least for polynomial $y$), the formula can be seen by linear algebra. One has part of the infinite dimensional matrix as follows: $$\left[\begin{array}{cccccccc}*&y(0)&y'(0)&y''(0)&y'''(0)&y^{(4)}(0)&y^{(5)}(0)&y^{(6)}(0)\\ A&1&\frac 1 2&\frac 1 6&\frac 1 {24}&\frac 1{120}&\frac 1{720}&\frac 1{5040}\\ B&0&1&\frac 1 2&\frac 1{6}&\frac 1{24}&\frac 1{120}&\frac 1{720}\\ C&0&0&1&\frac 1 2&\frac 1 6&\frac 1{24}&\frac 1{120}\\ D&0&0&0&1&\frac 1 2&\frac 1 6&\frac 1 {24}\\ E&0&0&0&0&1&\frac 1 2&\frac 1 6\\ F&0&0&0&0&0&1&\frac 1 2\\ G&0&0&0&0&0&0&1\end{array}\right]$$
By row operations to eliminate the higher-order components in $A$, one has $$y(0)=A-\frac 1 2B+\frac 1{12}C-\frac 1{720}E+\frac 1{30240}G,$$ modulo terms of order greater than $6$. This shows that $$y(1)=y(0)+(y(1)-y(0))=y(0)+B=A+\frac 1 2B+\frac 1{12}C-\frac 1 {720}E+\frac 1{30240}G+\cdots.$$ (So my comment has a small typo as well.)