Посмотрите на примере как обращаться с тем, что возвращает NDSolve[]
Код:
{Q1[s_], P1[s_], H1[s_]} = NDSolveValue[
{
q'[s] == p[s],
p'[s] == -q[s] - q[s]^2,
q[0] == 0.45,
p[0] == 0
},
{q[s], p[s], 1/2 p[s]^2 + 1/2 q[s]^2 + 1/3 q[s]^3},
{s, 0, 10},
Method -> {"SymplecticPartitionedRungeKutta",
"PositionVariables" -> {q[s]}}
] ;
{{Q2[s_], P2[s_], H2[s_]}} = {q[s], p[s],
1/2 p[s]^2 + 1/2 q[s]^2 + 1/3 q[s]^3} /. NDSolve[
{
q'[s] == p[s],
p'[s] == -q[s] - q[s]^2,
q[0] == 0.45,
p[0] == 0
},
{q[s], p[s], 1/2 p[s]^2 + 1/2 q[s]^2 + 1/3 q[s]^3},
{s, 0, 10},
Method -> {"SymplecticPartitionedRungeKutta",
"PositionVariables" -> {q[s]}} ] ;
ParametricPlot[{Q1[t], P1[t]}, {t, 0, 10}]
ParametricPlot[{Q2[t], P2[t]}, {t, 0, 10}]
Plot[{0, H1[t]}, {t, 0, 10}]
Plot[{0, H2[t]}, {t, 0, 10}]