There are have been quite a few questions on Mapleprimes, about numerically deriving additional information from the solutions returned by dsolve(..,numeric). That includes higher derivatives, integrals, inverting for particular function values, etc.
As it happens, I was having coffee last week with the principle developer of dsolve,numeric, and I asked him pointedly about all this. The answer I got back was pretty unequivocal: dsolve itself is the best thing to use in these situations, in general, because it will try to ensure that the requested accuracy (whether specified by default or explicit tolerance options) is attained for the requested quantities.
Consider fdiff (or evalf@D, which has a hook to fdiff). It will not compute as carefully as dsolve,numeric would. If dsolve,numeric detects difficulties near a requested point then it can automatically adjust stepsizes, etc. But that won't happen when using fdiff.
Let's compare the two methods given so far, pushed out only as far as the modest value of t=2.0. Fortunately we can compare against the exact result, which can be evalf'd at high working precision.
As we can see below, the result from dsolve,numeric itself is pretty much accurate up to its own default tolerances (which is incorporated into the procedures which it returns). The nested fdiff results, on the other hand, get wildly wrong as Digits changes.
solexact:=rhs(dsolve(sys union IC));
exp|- t |
Certainly there are easier ways, for lots of possible examples. But that is relatively unimportant, in terms of the general question.
Mathematical Software, Maplesoft