Easy way to plot ODE solutions from NDSolve?

Question

Inspired by the closed question Beautify a NDSolve Graph ! and a comment someone made to me not too long ago:

Is there some quick way to plot NDSolve results without going through the Plot and Evaluate[funcs /. sol] stuff?

Note the documentation for NDSolve is overflowing with examples of plotting solutions, via Plot and ParametericPlot, but perhaps there are other ways.

Examples

There is a variety of problems, but perhaps all can be addressed easily.

1. A simple ODE with a single solution:

var1 = {y};
ode1 = {y''[x] + y[x]^3 == Cos[x]};
ics1 = {y[0] == 0, y'[0] == 1};
sol1 = NDSolve[{ode1, ics1}, var1, {x, 0, 10}];    

2. A quadratic ODE with two solutions:

var2 = {y};
ode2 = {y''[x]^2 + y[x] y'[x] == 1};
ics2 = {y[0] == 0, y'[0] == 0};
sol2 = NDSolve[{ode2, ics2}, var2, {x, 0, 1}];

3. An ODE with a complex-valued solution:

var3 = {y};
ode3 = {y''[x] + (1 + Cos[x] I) y[x] == 0};
ics3 = {y[0] == 1, y'[0] == 0};
sol3 = First@NDSolve[{ode3, ics3}, var3, {x, 0, 20}];

4. A system of ODEs, with a single solution comprising multiple, real-valued, component functions:

var4 = {x1[t], x2[t], x3[t], x4[t]};
ode4 = {D[var4, t] == 
    Cos[Range[4] t] AdjacencyMatrix@
        CycleGraph[4, DirectedEdges -> True].var4 - var4 + 1};
ics4 = {(var4 /. t -> 0) == Range[4]};
sol4 = NDSolve[{ode4, ics4}, var4, {t, 0, 2}];

5. A vector-valued solution:

var5 = {x};
ode5 = {x'[t] ==
 (Cos[Range[4] t] AdjacencyMatrix@ CycleGraph[4, DirectedEdges -> True]).x[t] - x[t] + 1};
ics5 = {(x[t] /. t -> 0) == Range[4]};
sol5 = NDSolve[{ode5, ics5}, var5, {t, 0, 2}];

Answer 1

I wanted to share some undocumented techniques that give quick rough plots of NDSolve solutions. The keys points are this, the second one being quite handy at times:

• ListPlot[ifn] and ListLinePlot[ifn] wil plot an InterpolatingFunction ifn directly, if the domain and range are each real and one-dimensional. Points will be joined in line plots by straight segments; no recursive subdivision is performed.
• Similarly ListPlot[ifn'] and ListLinePlot[ifn'] will plot the derivatives of an InterpolatingFunction.
• The steps in the solution can be highlighted in line plots by either PlotMarkers -> Automatic or Mesh -> Full.
• One does not have to specify the domain for plotting, which is particularly useful when NDSolve runs into a singularity or stiffness and integration stops short. It's a convenient way to decide why the integration stopped.

The lack of recursive subdivision means ListLinePlot is good for examining the steps, but not good for examining the interpolation between the steps. The usual default interpolation order is 3, so the interpolation error is often a bit greater than the truncation error of NDSolve. For basic plotting, though, the steps by NDSolve are usually small enough that recursion is unnecessary to produce a good plot of the solution. If not, ListLinePlot[ifn, InterpolationOrder -> 3] will plot a smooth, interpolated path.

Normally, there's little difference between yfn = y /. First@NDSolve[..] and yfn = NDSolveValue[..], but see the second example. (For this reason and because the rules returned by NDSolve make it easy to substitute the solution into expressions such as invariants and residuals, I usually prefer NDSolve.)

Calls of the form NDSolve[..., y[x], {x, 0, 1}] result in ifn[x] instead of a pure InterpolatingFunction. To these, one has to apply Head to strip off the arguments in order to use ListPlot. See examples 3 and 5. (For this reason and because it difficult to substitute this form into y'[x], I usually prefer a call of the form NDSolve[..., y, {x, 0, 1}].)

Because ListLinePlot only plots real, scalar interpolating functions, complex-valued and vector-valued solutions are not as easily plotted as real, scalar interpolating functions. Some manipulation of the InterpolatingFunction is necessary. Perhaps someone else can come up with a better solution.

OP's examples:

1. Simple ODE

ListLinePlot[y /. First@sol1]

ListLinePlot[var1 /. First@sol1,  Mesh -> Full]
(* or ListLinePlot[y /. First@sol, PlotMarkers -> Automatic] *)

With the derivative:

ListLinePlot[{y, y'} /. First@sol1]

2. Nonlinear, multiple solutions

ListLinePlot[var2 /. sol2 // Flatten]
ListLinePlot[var2 /. #, PlotMarkers -> {Automatic, 5}] & /@ sol2 // 
 Multicolumn[#, 2] &
(* or ListLinePlot[y /. #, Mesh -> Full]& /@ sol // Multicolumn[#, 2]& *)

In this case, NDSolveValue is limited in what it does:

NDSolveValue[{ode2, ics2}, var2, {x, 0, 1}]

NDSolveValue::ndsvb: There are multiple solution branches for the equations, but NDSolveValue will return only one. Use NDSolve to get all of the solution branches.

3. Complex-valued solutions

This needs some extra handling so it is not as simple as applying ListLinePlot to the solution.

ListLinePlot[
 Transpose[{Flatten[y["Grid"] /. sol3], #}] & /@
  (ReIm[y["ValuesOnGrid"]] /. sol3), PlotLegends -> ReIm@y]

4. System with multiple components

If the call returned rules of the form x1 -> InterpolatingFunction[..] etc., mapping Head would not be needed. Otherwise, it would be simply pass a flat list of the interpolating functions. (The styling options are not really needed, of course.)

ListLinePlot[Head /@ Flatten[var4 /. sol4], PlotLegends -> var4, 
 PlotMarkers -> {Automatic, Tiny}, PlotStyle -> Thin]

5. Vector-valued solution

This, too, needs some extra manipulation of InterpolatingFunction.

ListLinePlot[
 Transpose[{Flatten[x["Grid"] /. sol5], #}] & /@
   (Transpose[x["ValuesOnGrid"]] /. First@sol5), 
 PlotLegends -> Array[Inactive[Part][x, #] &, 4]]

3D vector, with parametric plot:

var5b = x;
ode5b = {D[x[t], t] ==
  (Cos[Range[3] t] AdjacencyMatrix@ CycleGraph[3, DirectedEdges -> True]).x[t]};
ics5b = {x[0] == Range[-1, 1]};
sol5b = NDSolve[{ode5b, ics5b}, x, {t, 0, 2}];

ListPointPlot3D[x["ValuesOnGrid"] /. First@sol5b]
% /. Point[p_] :> {Thick, Line[p]}

Code dump

OPs code, in one spot, for cut & paste:

ClearAll[x,y,x1, x2, x3, x4];

(* simple ODE *)
var1 = {y};
ode1 = {y''[x] + y[x]^3 == Cos[x]};
ics1 = {y[0] == 0, y'[0] == 1};
sol1 = NDSolve[{ode1, ics1}, var1, {x, 0, 10}];

(* nonlinear, multiple solutions *)
ClearAll[y];
var2 = {y};
ode2 = {y''[x]^2 + y[x] y'[x] == 1};
ics2 = {y[0] == 0, y'[0] == 0};
sol2 = NDSolve[{ode2, ics2}, var2, {x, 0, 1}];

(* complex-valued solutions *)
var3 = {y};
ode3 = {y''[x] + (1 + Cos[x] I) y[x] == 0};
ics3 = {y[0] == 1, y'[0] == 0};
sol3 = First@NDSolve[{ode3, ics3}, var3, {x, 0, 20}];

(* system with multiple components *)
var4 = {x1[t], x2[t], x3[t], x4[t]};
ode4 = {D[var4, t] == 
    Cos[Range[4] t] AdjacencyMatrix@
        CycleGraph[4, DirectedEdges -> True].var4 - var4 + 1};
ics4 = {(var4 /. t -> 0) == Range[4]};
sol4 = NDSolve[{ode4, ics4}, var4, {t, 0, 2}];

(* vector-valued *)
ClearAll[x];
var5 = {x};
ode5 = {x'[
     t] == (Cos[Range[4] t] AdjacencyMatrix@
         CycleGraph[4, DirectedEdges -> True]).x[t] - x[t] + 1};
ics5 = {(x[t] /. t -> 0) == Range[4]};
sol5 = NDSolve[{ode5, ics5}, var5, {t, 0, 2}];

Answer 2

Whoever attempts an answer to this question has to repeat the details given in the mere documentation of NDSolve.

So what is important from this documentation page.

There are at first the different form that the NDSolve built-ins takes or requires the input.

There are three parts:

1. The equations of the ordinary differential problem.
2. The function u for the output or result.
3. The variable or set of variables and the ranges or domains of them.

The range, the domain and time dependence are therefore the main details of the possible plot of the resulting function u or set of resulting functions as the documentations page reports.

So the criteria of the first case are incomplete for a good specific answer.

1. A simple ODE with a single solution: But this is somehow close. This is because this tries to get close to the example "Second-order nonlinear ordinary differential equation:" from the documentation.

If to deal with an ordinary differential equations the start is the start and nowhere else. The CAS competitors has the macro

Identify the Type of an ODE

Then in the third step start the solution process.

1. Having the correct form of the ODE.
2. Analyse the type and correctness of the ODE.
3. Start the solution process.

From Maple there is this documentation page for the types of ODEs:

odeadvisor types

DEtools odeadvisor

To my first step

Description

•
Given an ODE, the odeadvisor command's main goal is to classify it according to standard text books (see dsolve,references), and to display a help page including related information for solving it (when the word help is given as an extra argument). The help pages include examples and Maple input lines, along with some advice, allowing you to adapt them to your problem. These help pages are also available by entering ?odeadvisor,; where is one of:

• First order ODEs

Abel,

Abel2A,

Abel2C,

Bernoulli,

Chini,

Clairaut,

dAlembert,

exact,

fully_exact_linear,

homogeneous,

homogeneousB,

homogeneousC,

homogeneousD,

homogeneousG,

linear,

patterns,

quadrature,

rational,

Riccati,

separable,

sym_implicit

In the page for "patterns" there is a discussion of the following ODE patterns:

y=g(y'),

x=g(y'),

0=G(x,y'),

0=G(y,y'),

y=G(x,y'),

x=G(y,y')

There is also a related parametric solving scheme.

Second order ODEs

Bessel,

Duffing,

ellipsoidal,

elliptic,

Emden,

erf,

exact_linear,

exact_nonlinear,

Gegenbauer,

Halm,

Hermite,

Jacobi,

Lagerstrom,

Laguerre,

Lienard,

Liouville,

linear_ODEs,

linear_sym,

missing,

Painleve,

quadrature,

reducible,

sym_Fx,

Titchmarsh,

Van_der_Pol

•
High order ODEs

quadrature,

missing,

exact_linear,

exact_nonlinear,

reducible,

linear_ODEs

After this typing two optimistic question of the solution strategy pathes are at first step solved: What is the potential effort for NDSolve, DSolve or simple integration with inverse functions? Is there a solution? Are the methods in need present in the built-ins of Mathematica. (I implicitly use the truth that Mathematica is better at ODEs than Maple, but not in all cases). If my ODE does not math one of the given types, there is chance I have to look up elsewhere.

This questions arise strict to the given effort of typing already and the regress to the types mentioned on the documentation page in the Mathematica documentation. The inherent strategy in the documentation of NDSolve is giving an exmaple to all solvable ODEs by NDSolves methods. So a well done paths to work with NDSolve is for a proper given ODE check the similarity with the examples in the documentation page of NDSolve.

For Mathemtica users there is a page already mentioned in other answers here that helps out:

solving ODE a work on published text book solutions

There is this answer on stackexchange:

how-to-determine-if-rhs-of-first-order-ode-is-separable-or-linear-or-neither-pa

Now starting with the very first example in the documentation page of Mathematica for NDSolve:

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}]

Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All]

Doing the work: Typing: This is an ordinary differential equation. Because it contain the function and the first order derivative of the function. It is first order. Only the function and the first derivative appear in the ode. This is a nonlinear ODE. This stems form the product of the function and the Cosine of the function. This is an nonlinear ODE with boundaries. There is a boundary conditions given: y[0]==1. There is a interval or domain. The interval is {x,0,30} in the Reals. Stemming from all the facts is that this nonlinear differential equation may be solved with NDSolve.

We did run the built-in with this special inputs.

Take a look at the solution function.

This is since this detail is given in the Mathematica documentation page for NDSolve:

"NDSolve gives results in terms of InterpolatingFunction objects. "

This is results are given in the representation of a rule or set of rules.

This detail implies there is not need for calculating another InterpolatingFunction with the results!

There are some historical inconviniencess with InterpolatingFunction: incompatible-interpolatingfunction-between-v9-and-earlier-and-v10.

So study the details of Your Mathematica version. InterpolatingFunction has some specifics that make it specially suitable for being the solution function for NDSolve.

• it works like Function.

• it interpolates and approximates.

• it incorporates the Domain, a superclass of Interval. Therefore it checks first whether the presented argument is in the domain and poses Interval arithmetics.

• one of the definitions of InterpolatingFunction in Wolfram Language is, that InterPolationFunction is the general results of NDSolve:

  NDSolve returns its results in terms of InterpolatingFunction objects.

It is possible to go deeper into details at this point of the discourse.

Our next step is to exammine the structure of the result of NDSolve.

To shorten we make use of:

s[[1, 1, 2]]

---

This part defintion gives us the InterpolatingFunction itself. With this we are able to work with the result of NDSolve as if is was an InterpolatingFunction object. Study the rest with another question.

s[[1, 1, 2]]["Methods"]

---

{"Coordinates", "DerivativeOrder", "Domain", "ElementMesh", \
"Evaluate", "GetPolynomial", "Grid", "InterpolationMethod", \
"InterpolationOrder", "MethodInformation", "Methods", \
"OutputDimensions", "Periodicity", "PlottableQ", "Properties", \
"QuantityUnits", "Unpack", "ValuesOnGrid"}

This is the list of Methods of the default InterpolatingFunction result of an arbitrary NDSolve input in Mathematica V12.0.0.

There is a closely related Methods regarding the given question: "PlottableQ":

s[[1, 1, 2]]["PlottableQ"]

---

True

This set of Methods makes up the information behind NDSolve work.

To decide whether to use a more basic or an advanced plotting built-in depends on the perspective. From the documentation page of InterpolatingFunction Plot is the built-in of recommendation.

This has many different advantages. One is nicely documented in the companion answer. Mathematica by default approximates functions values on grid for represenation purposes and interpolates for both function values in between given values of the x- and y-domain and the function curve itself. This methodologies are rather sophicticated for practical purposes. This not inside the scope of this answer. Mind Plot is the very best probably on the market at the moment of releases of new Mathematica in the new version. This is considered to be one of the major attractors for customers to buy the new version a priori without digging deep.

Another whow information is

s[[1, 1, 2]]["InterpolationMethod"]

---

"Hermite"

So the solution might differ from a closed representation. Hermite is method option presented in the documentation page for Interpolation.

"Interpolation supports a Method option. Possible settings include "Spline" for spline interpolation and "Hermite" for Hermite interpolation. "

While "Spline" is the already brute-force level, "Hermite" is more in and on the domain of the solution function set. Both rely on polynomials. This has the potential to be exact for polynomials but is ad hoc approximation for most mathematical functions. Polynomials span a vector space for themselves with infinite dimensionality and mightyness.

The mathematial truth for the given nonlinear ODE of first order with constant coefficients is, this is an exact ODE without separability. And it is not even of the type

x=G(y,y')

So there exists an integral for this problem.

The last equational identity exists only in points on the interval of definition in which the product y Cos[x+y] is unequal zero.

Solve[y Cos[x + y] == 0, {x, y}]

---

{{y -> 0}, {y -> -(\[Pi]/2) - x}, {y -> \[Pi]/2 - x}}

Reduce[y Cos[x + y] == 0, {x, y}]

---

(C1 [Element] Integers && (y == -([Pi]/2) - x + 2 [Pi] C1 || y == [Pi]/2 - x + 2 [Pi] C1)) || y == 0

So Reduce is better than Solve.

TrigExpand[ Cos[x + y]]

---

Cos[x] Cos[y] - Sin[x] Sin[y]

ties the sum of x and y together.

So we have an inifinite subset of the interval of existance for the solution in which there is no ad priori expectation for the existance of the solution. There is some filling for this gap allowing for the verification of the validity of the solution for the nonlinear and singular ODE. Both sides of the ODE have to be analysed for their properties carefully.

To shed some ideas on what selection Wolfram Inc. already made for us in the built-in NDSolve do this:

sa = AsymptoticDSolveValue[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, 
  y, {x, 0.00000000001, 10}]
Plot[{s[[1, 1, 2]][x], sa}, {x, 0, 3 \[Pi]/5}, 
 PlotRange -> {All, {0, 1.3}}]

This is some sort of hybrid between DSolve for exact solution and NDSolve for numerical solutions. This already works for a very long time.

The asymptotis solution is something really special:

1 + x Cos[1] + 1/2 x^2 (Cos[1]^2 - Sin[1] - Cos[1] Sin[1]) + 
 1/6 x^3 (-Cos[1] - 2 Cos[1]^2 - 3 Cos[1] Sin[1] - 4 Cos[1]^2 Sin[1] +
     Sin[1]^2 + Cos[1] Sin[1]^2) + 
 1/24 x^4 (-4 Cos[1]^2 - 11 Cos[1]^3 - 6 Cos[1]^4 + Sin[1] + 
    7 Cos[1] Sin[1] + 5 Cos[1]^2 Sin[1] - 6 Cos[1]^3 Sin[1] + 
    3 Sin[1]^2 + 13 Cos[1] Sin[1]^2 + 11 Cos[1]^2 Sin[1]^2 - 
    Sin[1]^3 - Cos[1] Sin[1]^3)-/+...

This shows that the sum in the cosine can be expanded and successive solutions exist. This is by far the more popular use in mathematics education for this example.

how-to-splice-together-several-instances-of-interpolatingfunction/19043#19043 offers some more details on the InterpolatingFunction capabilities not presented in the Mathematica documentation.

This Methods option Hermite still hides away a lot of the power of NDSolve for this problem.

Try for example this built-in combination:

Print /@ s[[1, 1, 2]] /@ {"Coordinates", "DerivativeOrder", "Domain", 
    "ElementMesh", Evaluate[], "Grid", "InterpolationMethod", 
    "InterpolationOrder", "Methods", "OutputDimensions", 
    "Periodicity", "PlottableQ", "Properties", "QuantityUnits", 
    "ValuesOnGrid"};
___
{{0.,0.00017069,0.00034138,0.00068276,0.00102414,0.00136552,0.00477932,0.00819312,0.0116069,0.0150207,0.0377779,0.0605351,0.0832922,0.106049,0.128807,0.168798,0.208789,0.238208,0.267626,0.297045,0.326464,0.355882,0.387732,0.419581,0.45143,0.483279,0.515129,0.546978,0.597245,0.647511,0.697777,0.748044,0.79831,0.848577,0.898843,0.967325,1.03581,1.10429,1.17277,1.24125,1.30974,1.37822,1.4467,1.54027,1.63384,1.7274,1.82097,1.91454,2.0081,2.10167,2.19524,2.2888,2.45347,2.61813,2.7458,2.87347,3.00114,3.1288,3.25647,3.38414,3.51181,3.63948,3.76714,3.89481,4.02248,4.15015,4.27782,4.40548,4.56912,4.73276,4.89639,5.02804,5.15969,5.29134,5.42299,5.55464,5.6641,5.77356,5.88301,5.99247,6.10193,6.12382,6.14571,6.1676,6.18949,6.21138,6.23327,6.25516,6.29895,6.34273,6.38651,6.4303,6.47408,6.51786,6.56164,6.62858,6.69552,6.76246,6.8294,6.89634,6.96328,7.03022,7.09715,7.16409,7.23103,7.29797,7.36491,7.43185,7.49879,7.56573,7.63267,7.72342,7.79857,7.85936,7.92015,7.98095,8.04174,8.10253,8.16332,8.22411,8.28491,8.37983,8.47475,8.56967,8.6646,8.75952,8.85444,8.94937,9.04429,9.16911,9.29392,9.41874,9.54355,9.66837,9.79319,9.918,10.0428,10.1865,10.3302,10.4739,10.6176,10.7614,10.9051,11.0488,11.1925,11.3362,11.5175,11.6988,11.8801,12.0614,12.2427,12.243,12.2432,12.2435,12.2438,12.244,12.2446,12.2451,12.2456,12.2509,12.2562,12.2615,12.2668,12.305,12.3432,12.3813,12.4195,12.4577,12.5216,12.5856,12.6495,12.7135,12.7774,12.8413,12.9053,12.9692,13.0331,13.0971,13.161,13.2249,13.2889,13.3528,13.4168,13.4898,13.5506,13.6114,13.6722,13.733,13.7938,13.8546,13.9284,14.0022,14.0599,14.1175,14.1752,14.2328,14.2905,14.3481,14.4058,14.4854,14.5651,14.6447,14.7243,14.804,14.8836,14.9632,15.0647,15.1663,15.2678,15.3693,15.4709,15.5724,15.6739,15.7755,15.9025,16.0295,16.1566,16.2836,16.4107,16.5377,16.6648,16.7918,16.9465,17.1012,17.256,17.4107,17.5654,17.7201,17.8749,18.0296,18.1843,18.339,18.4938,18.6485,18.7707,18.8929,18.9946,19.0963,19.1745,19.2528,19.331,19.4092,19.4874,19.5656,19.6438,19.722,19.8003,19.8785,19.9567,20.0349,20.1131,20.1913,20.2695,20.3478,20.426,20.5042,20.5824,20.6775,20.7726,20.8677,20.9628,21.0579,21.1531,21.2482,21.3433,21.4494,21.5555,21.6616,21.7677,21.8739,21.98,22.0861,22.1922,22.3129,22.4336,22.5542,22.6749,22.7956,22.9163,23.0369,23.1576,23.34,23.5225,23.6707,23.819,23.9673,24.1156,24.2639,24.4121,24.5604,24.7087,24.857,25.0052,25.1239,25.2426,25.3414,25.4403,25.5391,25.638,25.7368,25.8357,25.9345,26.0333,26.1322,26.231,26.3299,26.4287,26.5276,26.6264,26.7252,26.8241,26.9229,27.0218,27.1206,27.2195,27.3183,27.4172,27.516,27.6148,27.7137,27.8125,27.9238,28.0352,28.1465,28.2578,28.3691,28.4805,28.5918,28.7031,28.8529,29.0026,29.1524,29.3022,29.4519,29.6017,29.7515,29.8757,30.}}
___
0
___
{{0.,30.}}
____
None
___
{{0.},{0.00017069},{0.00034138},{0.00068276},{0.00102414},{0.00136552},{0.00477932},{0.00819312},{0.0116069},{0.0150207},{0.0377779},{0.0605351},{0.0832922},{0.106049},{0.128807},{0.168798},{0.208789},{0.238208},{0.267626},{0.297045},{0.326464},{0.355882},{0.387732},{0.419581},{0.45143},{0.483279},{0.515129},{0.546978},{0.597245},{0.647511},{0.697777},{0.748044},{0.79831},{0.848577},{0.898843},{0.967325},{1.03581},{1.10429},{1.17277},{1.24125},{1.30974},{1.37822},{1.4467},{1.54027},{1.63384},{1.7274},{1.82097},{1.91454},{2.0081},{2.10167},{2.19524},{2.2888},{2.45347},{2.61813},{2.7458},{2.87347},{3.00114},{3.1288},{3.25647},{3.38414},{3.51181},{3.63948},{3.76714},{3.89481},{4.02248},{4.15015},{4.27782},{4.40548},{4.56912},{4.73276},{4.89639},{5.02804},{5.15969},{5.29134},{5.42299},{5.55464},{5.6641},{5.77356},{5.88301},{5.99247},{6.10193},{6.12382},{6.14571},{6.1676},{6.18949},{6.21138},{6.23327},{6.25516},{6.29895},{6.34273},{6.38651},{6.4303},{6.47408},{6.51786},{6.56164},{6.62858},{6.69552},{6.76246},{6.8294},{6.89634},{6.96328},{7.03022},{7.09715},{7.16409},{7.23103},{7.29797},{7.36491},{7.43185},{7.49879},{7.56573},{7.63267},{7.72342},{7.79857},{7.85936},{7.92015},{7.98095},{8.04174},{8.10253},{8.16332},{8.22411},{8.28491},{8.37983},{8.47475},{8.56967},{8.6646},{8.75952},{8.85444},{8.94937},{9.04429},{9.16911},{9.29392},{9.41874},{9.54355},{9.66837},{9.79319},{9.918},{10.0428},{10.1865},{10.3302},{10.4739},{10.6176},{10.7614},{10.9051},{11.0488},{11.1925},{11.3362},{11.5175},{11.6988},{11.8801},{12.0614},{12.2427},{12.243},{12.2432},{12.2435},{12.2438},{12.244},{12.2446},{12.2451},{12.2456},{12.2509},{12.2562},{12.2615},{12.2668},{12.305},{12.3432},{12.3813},{12.4195},{12.4577},{12.5216},{12.5856},{12.6495},{12.7135},{12.7774},{12.8413},{12.9053},{12.9692},{13.0331},{13.0971},{13.161},{13.2249},{13.2889},{13.3528},{13.4168},{13.4898},{13.5506},{13.6114},{13.6722},{13.733},{13.7938},{13.8546},{13.9284},{14.0022},{14.0599},{14.1175},{14.1752},{14.2328},{14.2905},{14.3481},{14.4058},{14.4854},{14.5651},{14.6447},{14.7243},{14.804},{14.8836},{14.9632},{15.0647},{15.1663},{15.2678},{15.3693},{15.4709},{15.5724},{15.6739},{15.7755},{15.9025},{16.0295},{16.1566},{16.2836},{16.4107},{16.5377},{16.6648},{16.7918},{16.9465},{17.1012},{17.256},{17.4107},{17.5654},{17.7201},{17.8749},{18.0296},{18.1843},{18.339},{18.4938},{18.6485},{18.7707},{18.8929},{18.9946},{19.0963},{19.1745},{19.2528},{19.331},{19.4092},{19.4874},{19.5656},{19.6438},{19.722},{19.8003},{19.8785},{19.9567},{20.0349},{20.1131},{20.1913},{20.2695},{20.3478},{20.426},{20.5042},{20.5824},{20.6775},{20.7726},{20.8677},{20.9628},{21.0579},{21.1531},{21.2482},{21.3433},{21.4494},{21.5555},{21.6616},{21.7677},{21.8739},{21.98},{22.0861},{22.1922},{22.3129},{22.4336},{22.5542},{22.6749},{22.7956},{22.9163},{23.0369},{23.1576},{23.34},{23.5225},{23.6707},{23.819},{23.9673},{24.1156},{24.2639},{24.4121},{24.5604},{24.7087},{24.857},{25.0052},{25.1239},{25.2426},{25.3414},{25.4403},{25.5391},{25.638},{25.7368},{25.8357},{25.9345},{26.0333},{26.1322},{26.231},{26.3299},{26.4287},{26.5276},{26.6264},{26.7252},{26.8241},{26.9229},{27.0218},{27.1206},{27.2195},{27.3183},{27.4172},{27.516},{27.6148},{27.7137},{27.8125},{27.9238},{28.0352},{28.1465},{28.2578},{28.3691},{28.4805},{28.5918},{28.7031},{28.8529},{29.0026},{29.1524},{29.3022},{29.4519},{29.6017},{29.7515},{29.8757},{30.}}
___
Hermite
___
{3}
___
{Coordinates,DerivativeOrder,Domain,ElementMesh,Evaluate,GetPolynomial,Grid,InterpolationMethod,InterpolationOrder,MethodInformation,Methods,OutputDimensions,Periodicity,PlottableQ,Properties,QuantityUnits,Unpack,ValuesOnGrid}
___
{}
___
{False}
___
True
___
{Properties,Unpack}
___
{None,None}
___
{1.,1.00009,1.00018,1.00037,1.00055,1.00074,1.00257,1.00439,1.0062,1.008,1.01967,1.03078,1.04131,1.05123,1.06053,1.07533,1.08812,1.09622,1.10322,1.10912,1.11392,1.11765,1.12049,1.12212,1.12258,1.1219,1.12014,1.11732,1.11086,1.1021,1.09122,1.07841,1.06386,1.04775,1.03023,1.00441,0.976662,0.947338,0.916748,0.885177,0.852878,0.820082,0.786996,0.741656,0.69654,0.652011,0.608381,0.565922,0.524863,0.485398,0.447688,0.411861,0.353664,0.301905,0.266256,0.234456,0.206369,0.181804,0.160526,0.14227,0.126763,0.113728,0.102903,0.0940453,0.0869391,0.081401,0.0772818,0.0744679,0.0726511,0.0727831,0.0748944,0.0781192,0.0828351,0.0892124,0.0974674,0.107857,0.118326,0.130635,0.144942,0.16137,0.17998,0.183963,0.188032,0.192185,0.196421,0.200739,0.205136,0.209611,0.218785,0.228236,0.237938,0.247858,0.257957,0.268193,0.278518,0.294351,0.310059,0.325416,0.340186,0.35413,0.367018,0.378633,0.388785,0.397312,0.404091,0.409036,0.412102,0.413283,0.412607,0.410139,0.405966,0.39779,0.389038,0.380797,0.371638,0.361679,0.351034,0.339817,0.328139,0.316107,0.303819,0.284349,0.264816,0.245511,0.226677,0.208518,0.191196,0.174833,0.159516,0.141054,0.124538,0.109948,0.0972074,0.0862034,0.0768002,0.068851,0.0622087,0.0559982,0.0511424,0.0474779,0.0448741,0.0432353,0.042501,0.0426468,0.0436838,0.0456606,0.0496322,0.055498,0.063608,0.0743823,0.0882591,0.0882818,0.0883046,0.0883273,0.08835,0.0883727,0.0884182,0.0884638,0.0885093,0.0889664,0.0894265,0.0898896,0.0903557,0.093809,0.0974201,0.101191,0.105121,0.109212,0.116416,0.124055,0.132108,0.140546,0.149326,0.158394,0.167684,0.177115,0.186597,0.196026,0.205292,0.214277,0.22286,0.230921,0.238342,0.245904,0.251344,0.255936,0.259622,0.26236,0.264123,0.2649,0.264525,0.26274,0.260405,0.257288,0.253438,0.24891,0.243763,0.23806,0.231868,0.222634,0.212766,0.202435,0.191804,0.181023,0.17023,0.159547,0.146254,0.133495,0.121413,0.11011,0.0996517,0.0900721,0.0813788,0.0735577,0.0649513,0.0575749,0.0513254,0.0460933,0.0417705,0.038256,0.03546,0.0333063,0.0314627,0.0304007,0.0300739,0.0304702,0.0316118,0.0335564,0.0363973,0.0402634,0.0453162,0.0517421,0.0597363,0.0694748,0.0784794,0.0886511,0.0979611,0.107954,0.116023,0.124336,0.132792,0.14127,0.149633,0.157735,0.16542,0.172531,0.178917,0.184438,0.188974,0.192428,0.19473,0.19584,0.195748,0.194476,0.19207,0.188602,0.184164,0.177617,0.170007,0.161561,0.152504,0.143056,0.133419,0.123779,0.114295,0.104061,0.0943344,0.0852243,0.0768022,0.0691074,0.0621524,0.0559274,0.0504065,0.0449352,0.0402586,0.0363043,0.0330001,0.0302775,0.0280745,0.0263373,0.0250211,0.0237561,0.0232974,0.0234984,0.0242239,0.0255076,0.027406,0.0299991,0.0333889,0.037696,0.0430527,0.0495893,0.0574134,0.0646425,0.0727208,0.0800535,0.0878621,0.0960412,0.104448,0.112904,0.121198,0.129097,0.136356,0.142735,0.148015,0.15201,0.154585,0.155659,0.155213,0.153288,0.149977,0.145418,0.139783,0.133266,0.126074,0.11841,0.110471,0.102439,0.094472,0.0867066,0.079253,0.0713366,0.0640074,0.0573152,0.0512807,0.0459003,0.0411527,0.0370037,0.0334113,0.0293756,0.0261518,0.0236347,0.0217328,0.0203714,0.0194947,0.0190656,0.019036,0.0193018}

So plenty of information about our plot.

s[[1, 1, 2]]["MethodInformation"@#] &~Scan~s[[1, 1, 2]]["Methods"]

This give some internal representation overview and the complexity of the methods hidden in the InterpolatingFunction and in this case the result part of NDSolve.

This trick or insight from @carlwoll shows the flexibility of the combination I recommend:

ipf1 = Interpolation[Table[{x, Sin[x]}, {x, 0, 1, 0.1}]];
ipf2 = Interpolation[Table[{x, Sin[x]}, {x, 1, Pi, 0.1}]];

if = NDSolveValue[
    {y'[x] == Piecewise[{{ipf1'[x], x<1}}, ipf2'[x]], y[0] == ipf1[0]},
    y,
    {x, 0, 3.1}
];

if //OutputForm

---

InterpolatingFunction[{{0., 3.1}}, <>]

Plot[if[x], {x, 0, 3.1}]

This offers even more insight in our case into the internals of Methods option Hermite. This is a piecewise methodology on subintervals of the solution interval. So this example shows some very special power of NDSolve in a very nicely behaving case. That is very advertising. The answer Chaining extrapolation handlers from the referenced question: how-to-splice-together-several-instances-of-interpolatingfunction puts the hidden details into view.

The same concepts are very adaptable in ListPlot and ListLinePlot but there is need to put more feature in them for the higher degree of properness and exactness as shown in my answer.

So we learned NDSolve has the output InterpolatingFunction and this function is essentially better documentation in the page for Interpolation. In this there is the starting example an overlay with Show of Plot and ListPlot. The impression make the Plot by far more rugged an emphasises the subintervals of the interpolation. That is not that way for the results of NDSolve. So there is a problem for the use of ListPLot, ListLinePlot, ... for the results of NDSolve making that appear not as smooth as it really is. The internal methodes have to matched to each other.

The given examples:

1. A simple ODE with a single solution:

   var1 = {y}; ode1 = {y''[x] + y[x]^3 == Cos[x]}; ics1 = {y[0] == 0, y'[0] == 1}; sol1 = NDSolve[{ode1, ics1}, var1, {x, 0, 10}]

   Plot[Evaluate[y[x] /. sol1], {x, 0, 30}, PlotRange -> All]

Plot[sol1[[1, 1, 2]][x], {x, 0, 30}]

This sol1[[1, 1, 2]][x] is interpretable. 1 enters the first brace level of the solution function, the second 1 is the rule itself and the last 2 addresses the right hand side of the rule.

/. or ReplaceAll is shorter as input. The part methodology is closer to the newest innovation with knowledge representation of Wolfram Language. For the purposes of the question both are head-a-head speed. The part operator is closer to the methodological concepts of Plot stemming from the history of Mathematica. Both stay at the same level at which the ListPlot built-in were discarded. All leave it to Plot internals for speed and all the rest.

Back to the ODE. This is a nonlinear ordinary differential equation with inhomogeniety of second order and degree three in the polynomial nonlinearity. This is as a homogeneous ordinary differential equation an exact one.

var1 = {y};
ode1 = {y''[x] + y[x]^3 == 0};
ics1 = {y[0] == 0, y'[0] == 1};
sol1 = NDSolve[{ode1, ics1}, var1, {x, 0, 10}]
sol1e = DSolve[{ode1, ics1}, var1, {x, 0, 10}]
Plot[sol1[[1, 1, 2]][x], {x, 0, 10}]

{{y -> Function[{x}, 2^(1/4) JacobiSN[x/2^(1/4), -1]]}}

Plot[sol1[[1, 1, 2]][x], {x, 0, 10}]

Since this solution of NDSolve is the complete solution und this is an ordinary differential equation the homogenous and inhomogeneous solution are linear additive it is possible to display the pure inhomogenous part of the solution.

Plot[sol1[[1, 1, 2]][x] - sol1e[[1, 1, 2]][x], {x, 0, 10}]

By this the problem shall be solved.

2. A quadratic ODE with two solutions:

var2 = {y}; ode2 = {y''[x]^2 + y[x] y'[x] == 1}; ics2 = {y[0] == 0, y'[0] == 0}; sol2 = NDSolve[{ode2, ics2}, var2, {x, 0, 1}] Plot[{sol2[[1, 1, 2]][x], sol2[[2, 1, 2]][x]}, {x, 0, 1}]

3. An ODE with a complex-valued solution:

var3 = {y}; ode3 = {y''[x] + (1 + Cos[x] I) y[x] == 0}; ics3 = {y[0] == 1, y'[0] == 0}; sol3 = First@NDSolve[{ode3, ics3}, var3, {x, 0, 20}] ReImPlot[ImRe[sol3][[1, 1, 2]][x], {x, 0, 20}, PlotLegends -> Automatic]

4. A system of ODEs, with a single solution comprising multiple, real-valued, component functions:

var4 = {x1[t], x2[t], x3[t], x4[t]}; ode4 = {D[var4, t] == Cos[Range4 t] AdjacencyMatrix@ CycleGraph[4, DirectedEdges -> True].var4 - var4 + 1}; ics4 = {(var4 /. t -> 0) == Range4}; sol4 = NDSolve[{ode4, ics4}, var4, {t, 0, 2}];

Plot[{sol4[[1, 1, 2]], sol4[[1, 2, 2]], sol4[[1, 3, 2]], 
  sol4[[1, 4, 2]]}, {t, 0, 2}, PlotLegends -> Automatic]

A comparison addresses else:

Plot[{sol4[[1, 1, 2]], sol4[[1, 2, 2]], sol4[[1, 3, 2]], 
  sol4[[1, 4, 2]]}, {t, 0, 2}, PlotLegends -> var4, PlotStyle -> Thin]

5. A vector-valued solution:

var5 = {x}; ode5 = {x'[t] == (Cos[Range4 t] AdjacencyMatrix@ CycleGraph[4, DirectedEdges -> True]).x[t] - x[t] + 1}; ics5 = {(x[t] /. t -> 0) == Range4}; sol5 = NDSolve[{ode5, ics5}, var5, {t, 0, 2}];

Plot[sol5b[[1, 1, 2]][t], {t, 0, 2}]

ParametricPlot3D[sol5b[[1, 1, 2]][t], {t, 0, 2}]