For the following initial value problem

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


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:



  • $\begingroup$ @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. $\endgroup$
    – auxsvr
    Nov 26, 2012 at 17:41
  • $\begingroup$ eq[b_] := NDSolve[{x'[t] + x[t] == 0, x[0] == b}, x, {t, 0, 1}][[1]] and x'[1/2] /. eq[2] works? $\endgroup$
    – chris
    Nov 26, 2012 at 18:18
  • $\begingroup$ @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. $\endgroup$
    – m_goldberg
    Nov 26, 2012 at 18:30
  • $\begingroup$ @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? $\endgroup$
    – auxsvr
    Nov 26, 2012 at 20:33
  • $\begingroup$ 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] $\endgroup$
    – chris
    Nov 26, 2012 at 20:52

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.