I have played with the ODE code from Hull in 76 from Toronto -- but I am stuck -- I cannot see how to unwind one set of constructs
He has an infinite loop 9999
Before the loop he has a set of if statements controlled by an indicator flag ind , but the else if statements cross into the infinite loop. at ind == 3
he uses returns with the indicator flag to exit.
I was thinking the best way to pull it apart is a set of subroutines -- but someone else may have a bright idea.
Unfortunately they buried the RK method inside the main subroutine as well so it interesting
I enclose the code and the manual
John, there is a loop formed by statement 99999 and a GOTO to that label, but that loop is not infinite: there are six RETURN statements within that loop, any one of which serves to terminate the loop after it has been entered.
The results given by DVERK for your 2nd order ODE problem agree to at least four significant digits with that given by Ernst Hairer's DOPRI (4-5 pair) Dormand-Prince code at http://www.unige.ch/~hairer/software.html .
The DVERK code is old, and you may make it cosmetically better, if you wish, but there is no need to fix anything. Note that this code implements an RK code with two different integrators (5-6 pair) for step-size control based on error estimates, so it should be expected to be somewhat more complex than a simple RK-4 integrator with no step-size control. If there is anything "buried" in the code, it had to be. The integrator code is entirely within one subroutine; the only other subroutine that is needed to complete the code is a short subroutine that calculates the vector of derivatives of the dependent variables. This is the standard way of separating user-specified problem descriptions from the body of the integrator code.
Thank you for the comment. I agree the code is very good in terms of the results. It just looks old and messy, it does not need fixing, just making it more understandable to the average reader.