differential equations - Only final result from NDSolve

Finally, I started to play with differential equations in Mathematica.

And I have faced the problem, which seems to me so basic that I'm afraid this question is going to be closed soon. However, I've failed in finding solution here or in the documentation center.

My question is: how to set that NDSolve will not save whole InterpolationFunction for the result?

I'm only interested in final coordinates or for example each 100th. Is there simple way to achieve that?

---

Anticipating questions:

• 
  

  I know I can do something like r /@ Range[0.1, 1, .1] /. sol at the end, but still, there is whole interpolating function in memory. I want to avoid it because my final goal is to do N-Body simulation where N is huge and I will run out of the memory quite fast. What is important to me is only the set of coordinates as far in the future as it can be, not intermediate values.

  
  


• 
  

  I can write something using Do or Nest but I want to avoid it since NDSolve allows us to implement different solving methods in handy way.

  
  



• 
  

  I saw WolframDemonstrations/CollidingGalaxies and it seems there is an explicit code with Do :-/

  
  


• 
  

  Another idea would be to put NDSolve into the loop but this seems to be not efficient. Could it be even compilable?

  
  

---

Just in case someone wants to show something, here is some sample code to play with:

G = 4 Pi^2 // N ;


sol = NDSolve[{
               r''[t] == -G r[t]/Norm[r[t]]^3,
               r[0] == {1, 0, 0},
               r'[0] == {0, 2 Pi // N, 0}
              },
              r,
              {t, 0, 1}, Method -> "ExplicitRungeKutta", MaxStepSize -> (1/365 // N)
             ]


 ParametricPlot3D[Evaluate[r[t] /. sol], {t, 0, 1}]

  (* Earth orbiting the Sun. Units: Year/AstronomicalUnit/SunMass 
     in order to express it simply*)

Answer

Here is a solution inspired from tutorial/NDSolveStateData (Mathematica 8):

G = 4 Pi^2 // N;

stateData = 
 First[

    NDSolve`ProcessEquations[
        {  r''[t] == -G r[t]/Norm[r[t]]^3, 
           r[0] == {1, 0, 0},
           r'[0] == {0, 2 Pi // N, 0}
           },
        r,
        {t, 0, 1}, 
         Method -> "ExplicitRungeKutta",
         MaxStepSize -> (1/365 // N)]]


res = Table[
  NDSolve`Iterate[stateData, i/365];
  rAndDerivativesValues = NDSolve`ProcessSolutions[stateData, "Forward"];
  stateData = First @ NDSolve`Reinitialize[stateData, Equal @@@ Most[rAndDerivativesValues]];
  rAndDerivativesValues,
  {i, 1, 365, 7 (* 7 = each week*)}
   ] ;

rAndDerivativesValues is the value of r[t] and its derivative at each week.
Example :

Here is one revolution of the Earth around the Sun :

ListPointPlot3D[#[[1, 2]] & /@ res ]