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

Question

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.

Answer 1

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]}