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

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]

-
@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[0] == b}, x, {t, 0, 1}][[1]] and x'[1/2] /. eq[2] 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[0] == b}, x, {t, 0, 1}][[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][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 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][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 NDSolveReinitialize[]: $state = First[NDSolveReinitialize[$state, \[FormalX][0] == -3/2]];  from which you can propagate and evaluate once more: NDSolveIterate[$state, 5];
sol = NDSolveProcessSolutions[\$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 NDSolveReinitialize[]` whenever you need to change initial conditions.

-