Fehlberg - Lear Routines

[JohnNichols]

NASA in 1978 published a small report by Lear from TRW for the Shuttle Program on some New Nystrom Integrators. 

Nystrom is an extension to the standard Runge Kutta systems. 

Lear's work (attached) is quoted by a lot of the people who develop RK routines, but as one paper recently noted, there is not much available on a direct solution of the general second order equation. d2y/dx2 = f(t,x,dy/dx) equation (1)

Interestingly, Fybra in 1999 in an an update to a 1972 book on moving loads on bridges provided some samples of direct calculation of the general equation but only showed the resulting graphs. I was hoping to replicate Fybra's work to use on our bridges. It appears there is very limited implementations of equation (1). I started with ACM 670, which is from the Numerical Recipes library, but it only does dy/dx = f(t,x). I found an simple Fortran implementation of equation 1 from Jacob Williams (Degenerate Conic Web Site). It is a pretty simple implementation, but it was a good starting point. 

I now have it working in VS with the latest Intel Fortran, but I am challenged as to how to work out the error on each point.  The solution provides the correct answer for the problem solved by the ACM 670 code.  It appears that RK methods set the error as the difference in the solution between say 3 and 4 function integration - but that is a local error not necessarily a global error. Lear provides a three and four function sample for equation (1), which is possibly one way - computationally they are fast so that is not an issue. 

It is a long time since I studied RK, so I was wondering if there are any ideas this forum may have on determining errors? 

I enclose the code. 

Thanks

John