Here I have a system as follows:
$$\frac{dx}{dt}=a(by-x); \frac{dy}{dt}=rx-xz; \frac{dz}{dt}=(xy)^n-bz$$
Here $x, y$ and $z$ are positive real variables. All the parameters $a, r$ and $b$ are all positive real numbers and $n$ is a natural number.
How can I a make the phase portrait 3D of the system by varying the natural number $n$ from say $1$ to $100$?
3D phase portrait can be built by analogy with 2D using SliceVectorPlot3D[]
p[s1_, a1_, b1_, r1_, n1_] := Block[{s = s1, a = a1, b = b1, r = r1, n = n1}, v3D = {a*(b*y - x), r*x - x*z, (x*y)^n - b*z}; SliceVectorPlot3D[v3D/Norm[v3D], s, {x, -10, 10}, {y, -10, 10}, {z, -10, 10}, PlotTheme -> "Scientific", VectorColorFunction -> "BlueGreenYellow", VectorScale -> Small, VectorPoints -> Fine, PlotLabel -> Row[{"n=", n}], AxesLabel -> Automatic]] Table[p["XStackedPlanes", 2, 1, 4, n], {n, 1, 100, 33}] Table[p["YStackedPlanes", 2, 1, 4, n], {n, 1, 100, 33}] Table[p["ZStackedPlanes", 2, 1, 4, n], {n, 1, 100, 33}] 
a = 1; b = 1; r = 1; pf = ParametricNDSolveValue[{x'[t] == a*(b*y[t] - x[t]), y'[t] == r*x[t] - x[t]*z[t], z'[t] == (x[t]*y[t])^n - b*z[t], x[0] == x0, y[0] == x0, z[0] == x0}, {x[t], y[t], z[t]}, {t, 20}, {x0, n}]; Manipulate[ParametricPlot3D[pf[x0, n], {t, 0, 20}], {{x0, -1, "x0"}, -2, 2}, {{n, 1, "n"}, 1, 100}]