# Plot of phase portrait of system of 4D non linear differential equations

While numerically solving ths system of differential equations:

sol[w_, xc_, yc_, ac_, bc_] := NDSolve[{
x'[Z] == -(Log[10] (1/2 (3 x[Z]^3 + 9 x[Z]^2 b[Z] + x[Z] (-3 + a[Z]^2 + 9 b[Z]^2 - 9 y[Z]^2 (1 + 4 w b[Z]^2)) + b[Z] (-1 + a[Z]^2 + 3 b[Z]^2 - 3 y[Z]^2 (3 - 4 w + 12 w b[Z]^2))))),
y'[Z] == -(Log[10] (1/2 y[Z] (3 + 3 x[Z]^2 + a[Z]^2 + 6 (1 - 2 w) x[Z] b[Z] + 3 b[Z]^2 - 9 y[Z]^2 (1 + 4 w b[Z]^2)))),
a'[Z] == -(Log[10] (1/2 a[Z] (-1 + 3 x[Z]^2 + a[Z]^2 + 6 x[Z] b[Z] + 3 b[Z]^2 - 9 y[Z]^2 (1 + 4 w b[Z]^2)))),
b'[Z] == -(Log[10] (x[Z])),
x[7] == xc, y[7] == yc, a[7] == ac, b[7] == bc},
{x, y, a, b}, {Z, 7, -2}];


With initial conditions w=0.01, xc=3*10^-5, yc=6*10^-13, ac=Sqrt[0.999441], bc=38*10^-4

I've plotted the phase space trajectory but I'm struggling with plotting the phase portrait using StreamPlot. in order to see the flow of the vector field in the $$y-a$$ plane, with $$y$$ being in the x axis and $$a$$, in the y axis.

• I'm struggling to understand what this means: "in the $y-a$ plane, with $y$ being in the x axis and $a$, in the y axis". $y-a$ is a number (or scalar function), not a plane. Do you mean $y-a=0$? I don't know how $y$ can be in the $x$ axis ("on" the axis?), if $x$ and $y$ are independent variables. It seems unlikely that the 4D flow stays in a plane, although the flow could have an invariant subspace. And the 2D flow of the projection of the 4D vector field ≠ 2D projection of the 4D flow of the vector field. I'm unsure how you want to go from 4D to 2D. Commented May 2, 2023 at 17:11

We can use sections for different b,x as follows

sol[w_, xc_, yc_, ac_, bc_] :=
Module[{Z},
sols = NDSolve[{x'[
Z] == -(Log[
10] (1/2 (3 x[Z]^3 + 9 x[Z]^2 b[Z] +
x[Z] (-3 + a[Z]^2 + 9 b[Z]^2 -
9 y[Z]^2 (1 + 4 w b[Z]^2)) +
b[Z] (-1 + a[Z]^2 + 3 b[Z]^2 -
3 y[Z]^2 (3 - 4 w + 12 w b[Z]^2))))),
y'[Z] == -(Log[
10] (1/2 y[
Z] (3 + 3 x[Z]^2 + a[Z]^2 + 6 (1 - 2 w) x[Z] b[Z] +
3 b[Z]^2 - 9 y[Z]^2 (1 + 4 w b[Z]^2)))),
a'[Z] == -(Log[
10] (1/2 a[
Z] (-1 + 3 x[Z]^2 + a[Z]^2 + 6 x[Z] b[Z] + 3 b[Z]^2 -
9 y[Z]^2 (1 + 4 w b[Z]^2)))),
b'[Z] == -(Log[10] (x[Z])), x[7] == xc, y[7] == yc, a[7] == ac,
b[7] == bc}, {x, y, a, b}, {Z, -2, 7}]; sols[[1]]];

{w = 0.01, xc = 3*10^-5, yc = 6*10^-13, ac = Sqrt[0.999441],
bc = 38*10^-4};


Let generate data

s = sol[w, xc, yc, ac, bc];data = Table[{Z, b[Z], x[Z]} /. s, {Z, -2, 7, 1}];


To estimate range for y and a we use

Plot[Evaluate[{y[Z], a[Z]} /. s], {Z, -2, 7},
PlotLegends -> {"y", "a"}]


Finally we can evaluate

frame = Table[
StreamPlot[{-(Log[
10] (1/2 y[
Z] (3 + 3 x[Z]^2 + a[Z]^2 + 6 (1 - 2 w) x[Z] b[Z] +
3 b[Z]^2 - 9 y[Z]^2 (1 + 4 w b[Z]^2)))), -(Log[
10] (1/2 a[
Z] (-1 + 3 x[Z]^2 + a[Z]^2 + 6 x[Z] b[Z] + 3 b[Z]^2 -
9 y[Z]^2 (1 + 4 w b[Z]^2))))} /. {x[Z] -> data[[i, 3]],
b[Z] -> data[[i, 2]]}, {y[Z], 0, .6}, {a[Z], 0, 1},
FrameLabel -> {"y", "a"}, PlotLabel -> i], {i,
Length[data]}];


Note, that in this case we have small b,x, hence phase portrait in every section looks same, for example

{frame[[1]], Last[frame]}