Find the derivatives of arbitrary order of a function found using NDSolve

Question

For the following initial value problem

BC[b_]:=x[0]==b;
eq[b_]:=NDSolve[{equation==0, BC[b]}, x, {t,0,1000}];
body[t_,b_]:= x[t]/.eq[b];

one could define the derivative of body with regard to t as

Dbody[t_,b_]:=x'[t]/.eq[b];

However I'm trying to find a way to define this for derivatives of arbitrary order with regard to t, and make the definition of Dbody redundant. Simple invocations of ', D and Derivative failed; for example:

Derivative[1][body][1,1]

returns

body'[1,1]

Answer 1

One possibility is to use the StateData` framework of NDSolve[] for the purpose. I'll give a sketch on how one might go about it, and leave the encapsulation as a complete routine to you.

Start with an initial state:

$state = First @ NDSolve`ProcessEquations[{\[FormalX]'[t] + \[FormalX][t] == 0,
                                           \[FormalX][0] == 1}, \[FormalX], t];

(I have chosen the function name to be a formal symbol for safety.)

Whenever you need to propagate your solution, you only need to call NDSolve`Iterate[]:

NDSolve`Iterate[$state, 5]

After propagation (one can go either forward or backward, depending on the need), a call to NDSolve`ProcessSolutions[] yields an InterpolatingFunction[] object...

sol = NDSolve`ProcessSolutions[$state]
   {\[FormalX] -> InterpolatingFunction[{{0., 3.}}, <>]}

...which can now be used for evaluation:

{\[FormalX][2], \[FormalX]'''[3]} /. sol
   {0.135335, -0.0526585}

Why go through all that trouble of using state data, you ask? That's because changing initial conditions is a snap with this framework, by virtue of the function NDSolve`Reinitialize[]:

$state = First[NDSolve`Reinitialize[$state, \[FormalX][0] == -3/2]];

from which you can propagate and evaluate once more:

NDSolve`Iterate[$state, 5];
    sol = NDSolve`ProcessSolutions[$state];
{\[FormalX][2/3], \[FormalX]''[5/2]} /. sol
   {-0.770126, -0.122977}

The gist of it then, is to maintain a state of your DE, propagate if needed before evaluating, and use NDSolve`Reinitialize[] whenever you need to change initial conditions.