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

For the following initial value problem

BC[b_]:=x==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[body][1,1]


returns

body'[1,1]

• @m_goldberg: The title does not correspond to what I'm asking. I already get derivatives of arbitrary order wrt t if the initial conditions are constant. It is the case of two variables that I haven't been able to find. – auxsvr Nov 26 '12 at 17:41
• eq[b_] := NDSolve[{x'[t] + x[t] == 0, x == b}, x, {t, 0, 1}][] and x'[1/2] /. eq works? – chris Nov 26 '12 at 18:18
• @auxsvr. Your original title didn't correspond to your real question. I did my best to fix it. If you don't like what I did, re-edit it, but try to make it clearer this time. – m_goldberg Nov 26 '12 at 18:30
• @chris, it works, but it is not exactly what I want. Do you have any idea how to wrap this in a function of {t,b}, so that taking derivatives of it wrt t returns the correct result? – auxsvr Nov 26 '12 at 20:33
• does this do what you want ? f[t1_, c_] := Module[{eq}, eq[b_] := NDSolve[{x'[t] + x[t] == 0, x == b}, x, {t, 0, 1}][]; x'[t1] /. eq[c]]; f[1/2, 2] – chris Nov 26 '12 at 20:52

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.

$state = First @ NDSolveProcessEquations[{\[FormalX]'[t] + \[FormalX][t] == 0, \[FormalX] == 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 NDSolveIterate[]: NDSolveIterate[$state, 5]

After propagation (one can go either forward or backward, depending on the need), a call to NDSolveProcessSolutions[] yields an InterpolatingFunction[] object...
sol = NDSolveProcessSolutions[$state] {\[FormalX] -> InterpolatingFunction[{{0., 3.}}, <>]}  ...which can now be used for evaluation: {\[FormalX], \[FormalX]'''} /. 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 NDSolveReinitialize[]: $state = First[NDSolveReinitialize[$state, \[FormalX] == -3/2]];  from which you can propagate and evaluate once more: NDSolveIterate[$state, 5];

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