ParametricPlot with differential equations

Question

in the beginning of the code I have:

 θ[t_] := Array[Subscript[θ, #][t] &, n];

Then after defining the differential equations and the constrains in EQ I can solve them and Plot them like this:

 s = NDSolve[EQ, θ[t], {t, 0, 100}, Method -> {"EquationSimplification" -> "Residual"} ]

 Plot[Evaluate[θ[t] /. s], {t, 0, 10}, PlotRange -> All, PlotLegends -> Automatic]

But when I try a ParametricPlot, it no longer works:

 ParametricPlot[Evaluate[{θ[t] /. s, θ'[t] /. s]}, {t, 0, 10}, PlotRange -> All, PlotLegends -> Automatic]

Any idea why? (I think it is in the θ'[t], but I dont know why)

EDIT (all code)

n = 2;
g = 9.8;
m = {};
l = {};
EOM = Array[# &, n];

For[i = 1, i <= n, ++i,
   AppendTo[m, 0.1];
   AppendTo[l, 0.5];]

\[Theta][t_] := Array[Subscript[\[Theta], #][t] &, n];

Ec[t_] := 0.5*Sum[(n + 1 - i)*m[[i]]*(l[[i]]^2* (Subscript[\[Theta], i]'[t])^2 + Sum[2*l[[i]]*l[[k]]*Subscript[\[Theta], i]'[t] *Subscript[\[Theta], k]'[t] * Cos[Subscript[\[Theta], k][t] - Subscript[\[Theta], i][t]], {k, 1, i - 1}]), {i, 1, n}];

Ep[t_] := -g*Sum[  (n + 1 - i)*m[[i]]*l[[i]]*Cos[Subscript[\[Theta], i][t]]   , {i, 1, n}];

Lagrangiano[t_] = Ec[t] - Ep[t];

For[i = 1, i <= n, ++i,
EOM[[i]] = 
    D[D[Lagrangiano[t], Subscript[\[Theta], i]'[t]], t] - 
    D[Lagrangiano[t], Subscript[\[Theta], i][t]] == 0]


thetainicial = Array[# &, n];
Dthetainicial = Array[# &, n];
For[i = 1, i <= n, ++i,
   thetainicial[[i]] = Subscript[\[Theta], i][0] == Pi/4;
   Dthetainicial[[i]] = Subscript[\[Theta], i]'[0] == 0]
EQ = Join[EOM, thetainicial, Dthetainicial];

s = NDSolve[EQ, \[Theta][t], {t, 0, 100}, Method -> {"EquationSimplification" -> "Residual"} ]

Plot[Evaluate[\[Theta][t] /. s], {t, 0, 10}, PlotRange -> All, PlotLegends -> Automatic]

Everything fine until this:

Plot[Evaluate[\[Theta]'[t] /. s], {t, 0, 10}, PlotRange -> All, PlotLegends -> Automatic]

ParametricPlot[Evaluate[{\[Theta][t] /. s, \[Theta]'[t] /. s}], {t, 0, 10}, PlotRange -> All, PlotLegends -> Automatic]

Neither one works.

Answer 1

The following modifications provide what I think is the desired answer.

s = Flatten@NDSolve[EQ, θ[t], {t, 0, 100}, 
    Method -> {"EquationSimplification" -> "Residual"}];

to eliminate an extra {}. Then, as before, θ[t] is obtained by

Plot[Evaluate[θ[t] /. s], {t, 0, 10}, PlotLegends -> Automatic]

The other two plots require modified code:

Plot[Evaluate[D[θ[t] /. s, t]], {t, 0, 10}, PlotLegends -> Automatic] 
ParametricPlot[Evaluate[Transpose@{θ[t] /. s, D[θ[t] /. s, t]}], {t, 0, 10}, 
    AspectRatio -> 1, PlotLegends -> Automatic]

Transpose is needed to pair each θ'[t] with the corresponding θ[t],