These methods suggest the following questions:

When will they fail?
Can they be improved?
Are there better techniques available?

These methods obviously cannot be applied directly on an infinite interval, so if confronted by one you must do something to it, like changing the independent variable to make it into a finite interval problem before even thinking of using one.

They can and will fail as well for finite problems, in two general cases; if $f$ is such that its higher derivatives become large, you may have trouble attempting to use them.

Examples of such are $sin ⁡ ( 100000 x )$ and $sin ⁡ ( 1 x )$ for $x$ near 0. These are definite trouble for techniques like this.

In the former case, though $f$ is bounded, using small $N$ values will not give a reasonable picture of $f$ , and in the latter case, no $N$ value will do so sufficiently near 0.

There is another way you can run into a problem and that is when, though $f$ is well behaved, the solution you seek becomes infinite at some point in the interval you are concerned with.

That is not really as bad or strange a problem as you might expect. It does mess up these methods if you try to apply them naively though.

First it only really means that the reciprocal of $y$ , the function $1 y$ , goes through zero at the point at which $y$ becomes infinite. It is really no more surprising for $1 y$ to be zero than it is for it to be anything else, so this can easily happen with a non-linear differential equation.

Second and most important, if you know $y '$ you also know $( 1 y ) '$ , and can apply these methods to $1 y$ just as easily as you can to $y$ . If we make the definition: $z = 1 y ,$ we find, for $y ' = f ( x , y )$

$d z d x = d ( 1 / y ) d x = − z 2 d y d x = − z 2 f ( x , 1 / z )$

and we can apply the rules discussed to solve for $z ( b )$ instead of $y ( b )$ given $z ( a )$ .

In fact when $y$ is greater than 1, the derivative of $z$ will be smaller than that of $y$ and so you can expect to find it easier to approximate the changes in $z$ than in $y$ . When $y$ explodes and goes to infinity, $z$ wanders quietly near 0, and its behavior can easily be tracked by our methods.

Thus a cautious approach to solving equations where this can happen is to set up your favorite method for $y$ and for $z ( z = 1 y )$ simultaneously, and use the $y$ approach when the magnitude of $y$ is less than 1 and the $z$ one when the magnitude of $z$ is at most 1.

In our example, in which we have $y ' = x + y$ , the equation for $z$ becomes $z ' = − z 2 x − z$ , which is quite as easy to handle as the original equation.

In our computation we worked with the $y '$ equation, though $y$ was greater than 1 throughout the interval. We could probably have done even better by working with the $z$ equation.

Of course in our example, the original equation is linear in $y$ and the inhomogeneous term is linear in $x$ , so this problem will not arise.

Exercise 26.4 Try your favorite method on the differential equation $z ' = − z 2 x − z$ , with $z ( 0 ) = 1$ , and compare the reciprocal of the values you obtain for $z ( 1 )$ , that is, $1 z ( 1 )$ , with those you obtained earlier for $y$ .

There is much lore about numerical methods, and the subject is called numerical analysis.

It used to be a very dry field, because it is tedious and dull to learn a method that you can't possibly use, something like trying to read a cookbook when you have no kitchen to cook in.

Nowadays you can easily play with these methods, and they have dramatic application to real time control so that the field is now actually fun.
You will note that in every case, extrapolation using larger intervals only was able to improve the accuracy of estimates here, by factors on the order of a hundred thousand to a million, for each method.

Is extrapolation of the kind we have used the best we can do here to improve solutions?

These extrapolations have the great advantage that they are easy to perform, but they gain one added power of accuracy at the cost of doubling $N$ .
In fact, if you play your cards right, you can gain one added power of accuracy by increasing $N$ by 2, by choosing the best possible weights, for the points and using the correspondingly accurate rule for approximating f in each interval.
If you imagine $y$ written as a power series in $x$ , theoretically each new point could be used to eliminate the contribution from one higher power, so that with $N$ points it could be possible to produce a rule whose error would occur only from the $N$ -th and higher derivatives of $y$ .