Hello,
I need to solve the equations of motion (4 first order equations) for stars in the gravitational potential of the Galaxy. In order to do this I use the fifth order RK routines as given in NR. Thinks work fine when the equations are integrated up to a time of the order of 1Gyr.
Problems may occur when I try to integrate for a longer time, of the order of the age of the Galaxy (~15Gyr). In this case, the solution often seems not to converge using smaller and smaller values for the tolerance level epsilon; i.e. even changing epsilon of 15-20% (say, from 6e-9 to 5e-9) solutions may still change by a factor of 2 or similar. The same happens if I slightly change the trial stepsize h.
On the other hand, taking a too small epsilon produces a "stepsize underflow in rkqs" error message.
Have anyone some idea of where I am wrong? Thank you!
Cristiano

Without looking at the details, it's impossible to give a precise diagnosis. The generic therapy is as follows: The solutions may be running into cusps or other singularities or instability. Revert to graphical and local analytic methods. Try to devise a change of variables that eliminates the numerical problem.

Compute equilibrium points, if any, and Liapunov exponents. Local analytic solutions or perturbations theories. State-transition and parameter sensitivity matrices. But start with graphics to get an intuitive insight into the difficulties. One picture is worth 1000 words and maybe a terraflop or two.

Try to monitor constants of motion, i.e. energy of system.
It should be constant :) if not something is going bad.