The problem you're running into is that the usual InterpolatingFunction
solutions from NDSolve
have InterpolationOrder->3
. This means that 4 derivatives of your result will be identically 0.
Now, it is possible to raise the order of the solutions by choosing different methods, but this approach is limited by what methods are available to be used. Instead, I would augment the NDSolve
call so that it directly produces higher derivatives of your solution. For instance, given your initial NDSolve
example:
NDSolveValue[{y'[t] == y[t], y[0] == 1}, y[x], {t, 0, 1}]
we could augment it by one order to get:
sol = NDSolveValue[
{yp'[t] == yp[t], y'[t] == yp[t], y[0] == 1, yp[0] == 1},
{y, yp},
{t, 0, 2}
];
I used NDSolveValue
because I find it more convenient. Also, note that I had to give the derivative a different name. Here's a comparison of the two solutions:
Plot[{sol[[1]]'[x] - Exp[x], sol[[2]][x] - Exp[x]}, {x, 0, 2}, PlotRange->All]

The interpolating function corresponding to yp
gives a much more accurate result when compared to y'
. However, augmenting the NDSolve
/NDSolveValue
call is a bit tedious, so here's a function that does this in general:
raiseOrder[Inactive[nds:NDSolve|NDSolveValue][eqns_, v_, {x_, r__}, opts:OptionsPattern[]], order_] := Module[
{dpatt, eqn, ic, x0, max, low, neweqns},
dpatt = Replace[v, List[d__]:>Alternatives[d]];
(* Separate ODEs from initial/boundary conditions *)
{eqn, ic} = Lookup[{False, True}] @ GroupBy[eqns, FreeQ[x]];
(* Find initial/boundary points *)
x0 = Cases[ic, (dpatt|Derivative[_][dpatt])[z_] :> z, Infinity] //Union;
(* Boundary conditions are not supported *)
If[MatchQ[x0, {_}],
x0 = First @ x0,
Return[Inactive[NDSolve][eqns, v, {x, r}, opts]]
];
(* Find max order for each dependent variable v *)
max = GroupBy[
Cases[eqn, Derivative[n_][d:dpatt][x]:>{d, n}, Infinity],
First->Last,
Max
];
(* Differentiate ODEs requested number of times *)
neweqns = NestList[
D[#, x]&,
eqn,
order
];
(* Include needed additional first order equations. Also,
* include substitutes for the lower order equations in neweqns *)
low = KeyValueMap[
Table[derivative[n-1][#1]'[x] == derivative[n][#1][x], {n, #2-1}]&,
max + order
];
(* Add initial conditions related to the additional ODEs in neweqns *)
ic = Flatten @ Join[
Most @ neweqns,
ic
] /. Derivative->derivative /. x->x0;
(* Evaluate NDSolve/NDSolveValue *)
nds[
Flatten @ Join[
Take[neweqns, -1] /. Derivative[n_][y_]:>derivative[n-1][y]',
low,
ic
],
Flatten @ KeyValueMap[
Table[derivative[n][#1], {n, 0, #2+order-1}]&,
max
],
{x, r},
opts
] /. derivative -> Derivative
]
derivative[0][y_] := y
Here is your example computed to order 10:
sol = First @ raiseOrder[
Inactive[NDSolve][{y'[x]==y[x], y[0]==1}, y, {x, 0, 1}],
10
]

Finally, here is the power series expanded around 0:
Series[y[x], {x, 0, 10}] /. sol //TeXForm
$1.`+x+0.5` x^2+0.16666666666666666` x^3+0.041666666666666664` x^4+0.008333333333333333` x^5+0.001388888888888889` x^6+0.0001984126984126984`
x^7+0.0000248015873015873` x^8+\text{2.7557319223985893$\grave{ }$*${}^{\wedge}$-6} x^9+\text{2.755731922398589$\grave{ }$*${}^{\wedge}$-7}
x^{10}+O\left(x^{11}\right)$
Compare this to the actual series:
Series[Exp[1. x], {x, 0, 10}] //TeXForm
$1+x+0.5` x^2+0.16666666666666666` x^3+0.041666666666666664` x^4+0.008333333333333333` x^5+0.0013888888888888887` x^6+0.00019841269841269839`
x^7+0.000024801587301587298` x^8+\text{2.7557319223985884$\grave{ }$*${}^{\wedge}$-6} x^9+\text{2.7557319223985883$\grave{ }$*${}^{\wedge}$-7}
x^{10}+O\left(x^{11}\right)$
Of course, this isn't a very good example because the expansion around 0 doesn't really require solving an ODE. Instead, let's compute the power series around 1:
Series[y[x], {x, 1, 10}] /. sol //TeXForm
$2.718281861392969`+2.718281861392969` (x-1)+1.3591409306964846` (x-1)^2+0.45304697689882817` (x-1)^3+0.11326174422470704` (x-1)^4+0.022652348844941408`
(x-1)^5+0.0037753914741569016` (x-1)^6+0.0005393416391652716` (x-1)^7+0.00006741770489565895` (x-1)^8+\text{7.490856099517662$\grave{
}$*${}^{\wedge}$-6} (x-1)^9+\text{7.490856099517661$\grave{ }$*${}^{\wedge}$-7} (x-1)^{10}+O\left((x-1)^{11}\right)$
And the actual answer:
Series[Exp[x], {x, 1., 10}] //TeXForm
$2.718281828459045`+2.718281828459045` (x-1.`)+1.3591409142295225` (x-1.`)^2+0.45304697140984085` (x-1.`)^3+0.11326174285246021`
(x-1.`)^4+0.02265234857049204` (x-1.`)^5+0.0037753914284153404` (x-1.`)^6+0.0005393416326307629` (x-1.`)^7+0.00006741770407884536`
(x-1.`)^8+\text{7.490856008760597$\grave{ }$*${}^{\wedge}$-6} (x-1.`)^9+\text{7.490856008760595$\grave{ }$*${}^{\wedge}$-7}
(x-1.`)^{10}+O\left((x-1.`)^{11}\right)$
NDSolve[{y'[t] == y[t], y[0] == 1}, y'[t], {t, 0, 1}]
will solve for y'[t]. By solving for arbitrary derivatives, you might be able to use Taylor series expansion directly? $\endgroup$y''[t] /. First@NDSolve[{y'[t] == y[t], y[0] == 1}, y''[t], {t, 0, 1}] /. t -> 0
returns4.00037
, not one, as I would like it to... $\endgroup$y''[t] /. First@NDSolve[{y'[t] == y[t], y[0] == 1}, y''[t], {t, 0, 1}]
and thenPlot[%,{t,0,1}]
and the result looks correctish (I think your t->0 is being evaluated at the wrong place)? I agree with @Eckhard too: try FindSequenceFunction or something maybe. OK, I see what happens. The function is 4+ at t=0 (PlotRange->All helped) but it becomes correct even an infinitesimal distance from t=0. That's sort of reasonable since t=0 is a boundary. $\endgroup$Method->"ExplicitRungeKutta"
with a large enough"DifferenceOrder"
seems to do the trick. I might have to brush up on my numerical analysis skills to be better able to predict the effect of various methods and options. $\endgroup$