Consider solving the differential equation $y' = f(t,y)$ at time points $t_i$ with the respective solutions $y_i$. The following multi-step method is based on the integral relation:
$y(t_{i+1}) = y(t_{i-1}) + \int^{t_{i+1}}_{t_{i-1}}f(t,y)dt$
Derive an implicit multi-step method by approximating the integrand $f(t,y)$ with the polynomial interpolant using function values at $t_{i-2}, t_{i-1}, t_{i}, t_{i+1}$.
OK, so I know this multi-step method is called the Nystrom method. I also know that the method I come up with needs to have the following form:
$y_{i+1} = y_{i-1} + h(\alpha f_{i-2} + \beta f_{i-1} + \gamma f_{i} + \phi f_{i+1})$
Here $\alpha, \beta, \gamma$, and $\phi$ are constants I have to find. $h$ is just the timestep interval size, and $f_p = f(t_p,y_p)$ for some arbitrary $p$.
Beyond this, I'm completely stuck, and don't know how to even get started. Through only a few examples I found online, it seems what sometimes works is to use Taylor series to expand about the point $t_{i-1}$ and then somehow construct a system of equations for the constants to solve for them. Unfortunately the examples I've found haven't been really illustrative though. So when I tried this I got nowhere helpful pretty fast. I'd really appreciate any help I can get with this.
What you really need is to write down a quadrature rule for the desired integral:
Find $\alpha, \beta, \gamma, ...$ such that $$ \int_{t_{i-1}}^{t_{i+1}} f(t,y) dt \approx \alpha f_{i-2} + ... $$
To do that, you can use the interpolation polynomial $P(t,y)$ for $f(t,y)$ with respect to $t$ (in your case the degree of the polynomial will be 3 = number of points - 1). Then $$ \int_{t_{i-1}}^{t_{i+1}} f(t,y) dt \approx \int_{t_{i-1}}^{t_{i+1}} P(t,y) dt = \alpha f_{i-2} + ... $$ For interpolation polnomial you can use Lagrange form: https://en.wikipedia.org/wiki/Lagrange_polynomial