How can I plot the 3D diagram of x[t], x`[t] and a parameter ηb as the three axes?

Question

I have ploted a 2D diagram with ParametricPlot for the equations. The codes

Clear["Global`*"]
α = 110.; β = 55.; δ = 1.; μ1 = 18.; μ2 = 42.; μ = μ2/μ1;
deltap = .1;
inipoint = 2.;
tlength = 100.;
w[λ_, ξ_] := (-((μ1*α)/2) Log[1 - (λ^(-4) + 2*λ^2 - 3)/α]
  - (μ2*β)/2 Log[1 - (λ^-4*ξ^4 + 2 λ^2*ξ^-2 - 3)/β])/μ1;
dw[λ_, ξ_] = D[w[λ, ξ], λ];
ξin[λ_, ξ_, x_] = (1 + (λ^3 - 1) (x^3 - 1)^-1 (ξ^3 - 1))^(1/3);
f[λ_, ξ_, x_] = dw[λ, ξin[λ, ξ, x]]/(1 - λ^3);
sup[x_] := ((δ + x^3)/(1 + δ))^(1/3)

Get["NumericalDifferentialEquationAnalysis`"];

np = 11; points = weights = Table[Null, {np}];
intf[x0_, ξ0_] := Block[{y = x0, ξ1 = ξ0}, 
  Do[points[[i]] = 
    GaussianQuadratureWeights[np, y, sup[y]][[i, 1]], {i, 1, np}];
  Do[weights[[i]] = 
    GaussianQuadratureWeights[np, y, sup[y]][[i, 2]], {i, 1, np}];
  int = Sum[(f[λ, ξ1, y] /. λ -> points[[i]])*
     weights[[i]], {i, 1, np}]; int]

eq1 := x''[t] + (1/2 x'[t]^2 (3 - δ/x[t]^3 (1 + δ/x[t]^3)^(-4/3)
  - 3 (1 + δ/x[t]^3)^(-1/3)) + intf[x[t], ξ[t]]
  - deltap)/x[t]/(1 - (1 + δ/x[t]^3)^(-1/3)) == 0;
eq2 := (3 ηb*(1 - (x[t]^-4*ξ[t]^4 + 2 x[t]^2*ξ[t]^-2 - 3)/β)) ξ'[t] ==
  ξ[t]*(μ (x[t]^2*ξ[t]^-2 - x[t]^-4*ξ[t]^4));

pfun = ParametricNDSolveValue[{eq1, eq2,
  ξ[0] == 1, x'[0] == 0, x[0] == inipoint},
  {x[t], x'[t]}, {t, 0, tlength}, {ηb}];

ParametricPlot[pfun[0.1], {t, 0, tlength}]

the result is

Then, I want to add the third axis (ηb from 0.1 to 1.) in 3D diagram, any ideas?

Answer 1

Using ηb as the vertical axis value

ParametricPlot3D[Evaluate[Table[{pfun[p][[1]], pfun[p][[2]], p}, {p, Range[.1, 1., .05]}]], 
 {t, 0, tlength}, BoxRatios -> 1, 
 PlotPoints -> 3000, MaxRecursion -> 10, 
 PlotLegends -> (Row[{ "p = ", #}] & /@ Range[.1, 1., .05]), 
 PlotRange -> All]