An Unstable Method of Solving Differential Equations

Here we are going to develop a method of solving differential equations which at the outset looks very good, but actually is very bad.  Let us develop an order 4 method to solve the equation
        (dy/dt)=f(t,y)
To do this, we will use an approximate formula for y(t[i]) from Chapter 4.1.  We don't calculate it, but the error is O(h^4).

list={{ti+h,y[ti+h]},{ti,y[ti]},{ti-h,y[ti-h]},{ti-2h,y[ti-2h]},{ti-3h,y[ti-3h]}};
p[t_]=Expand[InterpolatingPolynomial[list,t]];
Together[p'[ti]]

So we can substutite ti+h,ti, etc with w[i+1], w[i], etc, and we get the method:
        ((3(Subscript[w, i+1])+10(Subscript[w, i])-18(Subscript[w, i-1])+6(Subscript[w, i-2])-(Subscript[w, i-3]))/(12h))=f((Subscript[t, i]),(Subscript[w, i]))
that is
        (Subscript[w, i+1])=-(10/3)(Subscript[w, i])+6(Subscript[w, i-1])-2(Subscript[w, i-2])+(1/3)(Subscript[w, i-3])+4h f((Subscript[t, i]),(Subscript[w, i]))

It is not so hard to see that this is a method of order 4, but let us not dwell on that.  To illustrate how poor this method is, we will consider the differential equation
        (dy/dt)=0
        y(0)=1
So (Subscript[w, 0])=1.  To apply this method we need to also calculate (Subscript[w, 1]), (Subscript[w, 2]), and (Subscript[w, 3]).  We might use say a Runge-Kutta method to do this.  We will suppose that a little bit of numerical error sets in.  We will set h=0.0001.

w[0]=1;
w[1]=1;
w[2]=1;
w[3]=1.0000000000001;
Do[w[i+1]=-10/3w[i]+6w[i-1]-2w[i-2]+1/3w[i-3],{i,3,30}]
TableForm[Table[{0.0001i,w[i]},{i,0,30}]]

As we can see, the answer becomes useless at about the 22nd step.  This is irespective of the value of h.  Why is this?  Well one way to solve this recurrence relation is by considering the roots of the polynomial:

NSolve[s^4==-10/3s^3+6s^2-2s+1/3,s]

So the solution to the recurence relation will include a term of the form A(-4.7028)^i where A is some very small constant.  Ideally we should have A=0, but small numerical errors will obtain that A is a small number.  But no matter how small A is, after a few steps, this term will begin to dominate the answer.

Hence the upshot of all this is the the numerical method we have developed is unstable, and as such, is completely useless.