How to solve a linear difference equation? 
We have the linear multi-step process
$-y_k-y_{k+1}+y_{k+2}+y_{k+3}=4hf(t_k,y_k)$
for the initial value problem $y'(t)=f(t,y(t)), y(0)=0$
Give every solution of the homogeneous  difference equation.

[Note, that this is just one part of the task, I am asking for help in solving the difference equation]
So I want to solve the homogeneous  difference equation:
$-y_k-y_{k+1}+y_{k+2}+y_{k+3}=0$
I have never done this, so I am not exactly sure on what to do, and I might need some help.
As far as I know, I have to convert this into the equation:
$-1-x+x^2+x^3=0$ and solve for $x$.
This is easy:
It is $x^3+x^2-x-1=(x-1)(x+1)^2$ and we have the roots $x=1$ and $x=-1$ (which counts double)
What to do now?
Do I have to find polynomials $p_1(x)$ and $p_2(x)$ such that:
$x^3+x^2-x-1=p_1(x)\cdot 1+p_2(x)\cdot (-1)^2$
This would give many solutions.
Thanks in advance for your help.
 A: The solution of the general third order difference equation is
$$
y_k=c_1q_1^k+c_2q_2^k+c_3q_3^k
$$
and for the case of a double root $q_2=q_3$
$$
y_k=c_1q_1^k+c_2q_2^k+c_3kq_2^k.
$$

In the concrete case you thus get
$$
y_k=c_1+(-1)^k(c_2+c_3k).
$$

In each step of the ODE integration there is an error contribution of truncation error and floating point rounding noise. This gets distributed over all the coefficients, so that the accumulated error grows from this alone. Now note that the factor $k$ is the number of steps, that is $t/h$, so that the third component grows inverse to the step size when comparing numerical approximations of different step sizes.

Numerical experiments for the method
Apply the method to the slightly non-linear problem
$$
y'(t)=p'(t)-0.1(\sin(0.01y(t))-\sin(0.01p(t))
$$
with the exact solution $y(t)=p(t)=\cos(t)$. For the initialization of the method, use the exact values $y_k=p(t_k)$ for $k=0,1,2$. Then plot solutions and their errors, scaled by $h$ as this is a first order method. This gives the graphs

One can see from the reduction of the error in the first plot resp. the uniform boundedness of the error in the second plot that this is indeed a first-order method. However, the rapid growth of the error for even this mild equation makes this rather unsuitable for practical applications.
