Determine the variables so as to make differential equation has global error of order 3 
Determine $\alpha, \beta$ and $\gamma$ so that the linear, multistep
method
$$ y_{j+4} - y_j + \alpha (y_{j+3} - y_{j+1}) = h [ \beta (f_{j+3} -
 f_{j+1}) + \gamma f_{j+2} ] $$
for the d.e. $y' = f(x)$ has global error of order $3$. Is the
resulting method stable?

Attempt
We assume $f_{j} = f(x_j,y_j)$. At the n+1 step, the solution (exact) is $y(t_n)$ and the approximation is $y_{n+1}$. So the way we find truncation error is by using taylor expansion on $y(t_n) - y_{n+1}$ and this works fine for single-step method. But, in this situation, how would we do it?
 A: The symmetry of the question suggests centering everything at $x_{j+2}$.  We can translate the question by making $x_{j+2}=0$ to ease the typing.  Note that $y_1=y(-h)$ and so on.  First it must be exact for $f(x)=1,y(0)=0,y=x$, so plug that in.  We get
$$4h+2\alpha h=h\gamma\\ \gamma=4+2\alpha$$
Then it must be exact for $f(x)=x,y(x)=\frac 12x^2$, so we have
$$0=h\beta(2h)\\\beta=0$$
It must be exact for $f(x)=x^2,y(x)=\frac 13x^3$, so we have
$$\frac {16}3h^3+\alpha\frac 23h^3=0\\ \alpha=-8\\ \gamma=-12$$
We are only looking at the odd part of $y$ around $x_{j+2}$ and the odd part of $x$.  It probably satisfies your definition of stability as the right side is zero, but it does not solve the equation.
A: Another test problem to get more coefficient identities at once it the classical $f(t,y)=y$ with solution $y(t)=Ce^x$. Select the constant so that $y_{j+2}=1$, then the equation reads
$$
2\sinh(2h)+2α\sinh(h)=h(2β\sinh(h)+γ)+O(h^{p+1})
$$
Now the left side is odd in $h$, so the right side can not have even terms, requiring $β=0$.
Inserting the power series of the exponential and comparing coefficients of equal degree gives next
\begin{align}
\left[2h+\frac{(2h)^3}{6}+\frac{(2h)^5}{5!}+...\right]+α\left[h+\frac{h^3}{6}+\frac{h^5}{5!}+...\right]&=\frac{hγ}2+O(h^{p+1})
\\[1em]\hline
2+α&=\frac{γ}2\\
\frac{8+α}6&=0
\end{align}
so that $α=-8$, $γ=-12$ and $p=4$.

As for stability, look at the characteristic polynomial for the left side,
$$
z^4-8z^3+8z-1=(z^2-1)(z^2-8z+1)=(z-1)(z+1)(z-4-\sqrt{15})(z-4+\sqrt{15})
$$
The root $z=-1$ gives a weakly unstable component, the root $z=4+\sqrt{15}$ makes the method strongly unstable, errors due to the initialization of the multi-step method and floating point truncations will be rapidly magnified.

possible typo in given equation?
Assuming the first component of the right side should lead to a non-trivial contribution, assume that there is a typo and it is a sum, not a difference. Then the test equation leads to
$$
2\sinh(2h)+2α\sinh(h)=h(2β\cosh(h)+γ)+O(h^{p+1})
$$
with the power series expansion
\begin{align}
\sum_k(2^{2k+1}+α)\frac{h^{2k+1}}{(2k+1)!}
&=hβ\sum_k\frac{h^{2k}}{(2k)!}+\frac{hγ}2+O(h^{p+1})
\\[1em]\hline
2+α&=β+\frac{γ}2\\
\frac{8+α}6&=\frac{β}2\\
\frac{32+α}{120}&=\frac{β}{24}
\end{align}
giving $β=12$, $α=28$, $γ=36$ with $p=6$ but even more instability.
