Export NDSolve data

Question

I have calculated a solution of a system of differential equations using NDSolve and generated a plot. How could I extract the data of the solutions? The input is as follows:

κ := 2; n := 1; m := 1 ; γ := 0.1 ;

NDSolve[
{Subscript[U, 2]''[t] + (n^2 + 5 m^2/Cos[m t]^2) Subscript[U, 2][t] == 0, 
 Subscript[F, 2]''[t] + (n^2 + 5 m^2/Cos[m t]^2) Subscript[F, 2][t] + 
     24 γ κ^3 m Sin[m t] Cos[m t]^2 Subscript[F, 8][t] ==
     0, 
Subscript[F, 8]''[t] + (n^2 + m^2) Subscript[F, 8][t] + 24 γ κ^3 m Sin[m t] Cos[m t]^2 Subscript[F, 2][t] == 0, 
Subscript[F, 2][0] == 0, Subscript[F, 2]'[0] == 1, Subscript[F, 8][0] == 0, Subscript[F, 8]'[0] == 1, Subscript[U, 2]'[0.] == 1, Subscript[U, 2][0] == 0}, {Subscript[U, 2], Subscript[F, 2], Subscript[F, 8]}, {t, 0, 0.4999 π}];

Plot[Evaluate[{Subscript[F, 2][t], Subscript[U, 2][t]} /. %], {t, 0, 0.4999 π}]

Answer 1

As noted by J. M. in the comments, the easiest way to do this is to create a Table of x and t values. However, if you want to actually extract the data points calculated by Mathematica, the InterpolatingFunctionAnatomy package allows you to get this data out. An example:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
soln = NDSolve[{x''[t] == -5 (1 + 1000 Exp[-100 t^2]) x[t], x[-5] == 0, x'[-5] == 1}, x, {t, -5, 5}]
xsoln = x /. First[First[soln]]
tvals = First[InterpolatingFunctionCoordinates[xsoln]];
xvals = InterpolatingFunctionValuesOnGrid[xsoln];
ListPlot[Transpose[{tvals, xvals}], PlotRange -> All]

Note the differing spacing of the data points in the quickly-varying region, due to Mathematica's adaptive step-size algorithms.

Answer 2

First of all, the kind of interpolation produced by NDSolve for an ODE can depend on the Method and the setting of InterpolationOrder. I will answer just for the OP's case, which is the usual case. NDSolve stores more data than just the abscissae and ordinates. It also stores derivative values at each abcissa, up to the order of the ODE. I will show how to package this in a table of the form $$\{ \,\{\,x_i,y_i,y_i',y_i'',\dots\,\},\dots\,\}$$

The easy way to do this is the following if sol is the solution to the OP's NDSolve command:

With[{f = Subscript[U, 2] /. First@sol},
 data = Table[{x, f[x], f'[x], f''[x]}, {x, First@f["Coordinates"]}]]

(*  {{0., 0., 1., -2.1684*10^-19}, {0.000107875, 0.000107875, 1., -0.00064725},...  *)

There are sometimes very slight round-off errors compared with the data stored in the solution. The following extracts the data from the solution's InterpolatingFunction for Subscript[U, 2]. Some description of the internal structure may be found here, What's inside InterpolatingFunction[{{1., 4.}}, <>]?, and some more can be found here, Interpolating data with a step. The usual InterpolatingFunction produced by NDSolve has the data stored in an argument of the form

{Developer`PackedArrayForm, idcs_, vals_}

where vals is a flattened array of the ordinates and derivative values. The indices idcs indicate where to partition vals to recover the array. Usually in a solution to an ODE these are equally spaced. The code below assumes all this. If it fails, then the solution at hand is not of the form I have been calling "usual."

solvals = 
  First@ Cases[
    Subscript[U, 2] /. First@sol,
    {Developer`PackedArrayForm, idcs_, vals_} :> 
     Partition[vals, Differences@idcs /. {n_Integer ..} :> n], 
    Infinity];
data = Transpose@ Join[Subscript[U, 2]["Coordinates"] /. First@sol, Transpose[solvals]]

(*  {{0., 0., 1., 0.}, {0.000107875, 0.000107875, 1., -0.00064725},...  *)