A high order method for the numerical solution of two-point boundary value problems

A one step finite difference scheme of order 4 for the numerical solution of the general two-point boundary value problemy=f(t,y),a t b, withg(y(a),y(b))=0 is presented. The global discretization error of the scheme is shown, in sufficiently smooth cases, to have an asymptotic expansion containing even powers of the mesh size only. This justifies the use of Richardson extrapolation (or deferred correction) to obtain high orders of accuracy. A theoretical examination of the new scheme for large systems of equations shows that for a given mesh size it generally requires about twice as much work as the Keller box scheme. However, the expectation of higher accuracy usually justifies this extra computational effort. Some numerical results are given which confirm these expectations and show that the new scheme can be generally competitive with the box scheme.