Solution (source code)

= Solution

Insert a smooth exact solution, so $f(y(t))=y'(t)$ and $g(y(t))=y''(t)$. Expand about $t_n$. The residual, left side minus right side, has the formal differential-operator symbol
$$
\mathcal R(z)=e^z-e^{-z}-\frac7{15}z(e^z+e^{-z})-\frac{16}{15}z+\frac1{15}z^2(e^z-e^{-z}),\qquad z=hD.
$$
Its odd symmetry makes every even coefficient vanish. The coefficients of $z,z^3,z^5$ are respectively
$$
2-\frac{14}{15}-\frac{16}{15}=0,\qquad
\frac13-\frac7{15}+\frac2{15}=0,\qquad
\frac1{60}-\frac7{180}+\frac1{45}=0.
$$
The next coefficient is $1/4725\ne0$, giving
$$
\mathcal R(hD)y(t_n)=\frac{h^7}{4725}y^{(7)}(t_n)+O(h^9).
$$
Hence the <local truncation error> has exact defect order seven and \b[the method has order six]. The zero-step <characteristic polynomial> is $\zeta^2-1$, whose two unit roots are simple, so it is also <zero-stable>. With sufficiently accurate starting values and exact $g=f'f$ evaluations, the usual accumulation of defects gives sixth-order <global error>. This is the <symmetric sixth-order two-derivative method>; approximate evaluations of $g$ must preserve the required accuracy if that order is to be retained.