ParametricNDSolve for a system of 4 non-linear ODEs and plotting the solutions

Question

I have four coupled first-order non-linear differential equations, denoted as: A1[x], A2[x], A3[x], A4[x] which are all functions of x. I have the following code which attempts to solve the equations using ParametricSolveND by varying one of the initial conditions of the parameter (namely, A4[0] which I have denoted as the parameter j).

ω1 = 2 π*5*10^9;
ω2 = 2 π*5*10^9;
ω3 = 2 π*3*10^9;
ω4 = ω1 + ω2 - ω3;
Cj = 329*10^-15;
LL = 100*10^-12;
a = 10*10^-6;
I0 = 3.29*10^-6;
CC0 = 39*10^-15;
k1 = (Sqrt[CC0 LL] *(ω1))/(a Sqrt[1 - Cj LL *(ω1)^2]);
k2 = (Sqrt[CC0 LL] *(ω2))/(a Sqrt[1 - Cj LL *(ω2)^2]);
k3 = (Sqrt[CC0 LL] *(ω3))/(a Sqrt[1 - Cj LL *(ω3)^2]);
k4 = (Sqrt[CC0 LL] *(ω4))/(a Sqrt[1 - Cj LL *(ω4)^2]);
Δkl = k1 + k2 - k3 - k4;
κ1 = (a^4*k1*k2*k3*k4*(k3 + k4 - k2))/(8*CC0*I0^2*LL^3*ω1^2);
κ2 = (a^4*k1*k2*k3*k4*(k3 + k4 - k1))/(8*CC0*I0^2*LL^3*ω2^2);
κ3 = (a^4*k1*k2*k3*k4*(k1 + k2 - k4))/(8*CC0*I0^2*LL^3*ω3^2);
κ4 = (a^4*k1*k2*k3*k4*(k1 + k2 - k3))/(8*CC0*I0^2*LL^3*ω4^2);
α11 = (a^4*k1^5)/(16*CC0*I0^2*LL^3*ω1^2);
α12 = (a^4*k1^3*k2^2)/(8*CC0*I0^2*LL^3*ω1^2);
α13 = (a^4*k1^3*k3^2)/(8*CC0*I0^2*LL^3*ω1^2);
α14 = (a^4*k1^3*k4^2)/(8*CC0*I0^2*LL^3*ω1^2);
α21 = (a^4*k2^3*k1^2)/(8*CC0*I0^2*LL^3*ω2^2);
α22 = (a^4*k2^5)/(16*CC0*I0^2*LL^3*ω2^2);
α23 = (a^4*k2^3*k3^2)/(8*CC0*I0^2*LL^3*ω2^2);
α24 = (a^4*k2^3*k4^2)/(8*CC0*I0^2*LL^3*ω2^2);
α31 = (a^4*k3^3*k1^2)/(8*CC0*I0^2*LL^3*ω3^2);
α32 = (a^4*k3^3*k2^2)/(8*CC0*I0^2*LL^3*ω3^2);
α33 = (a^4*k3^5)/(16*CC0*I0^2*LL^3*ω3^2);
α34 = (a^4*k3^3*k4^2)/(8*CC0*I0^2*LL^3*ω3^2);
α41 = (a^4*k4^3*k1^2)/(8*CC0*I0^2*LL^3*ω4^2);
α42 = (a^4*k4^3*k2^2)/(8*CC0*I0^2*LL^3*ω4^2);
α43 = (a^4*k4^3*k3^2)/(8*CC0*I0^2*LL^3*ω4^2);
α44 = (a^4*k4^5)/(16*CC0*I0^2*LL^3*ω4^2) // N;

system = {A1'[x] == I*κ1*Conjugate[A2[x]]*A3[x]*A4[x]*E^(-I*Δkl*x) + I*A1[x]*(α11*Abs[A1[x]]^2 + α12*Abs[A2[x]]^2 + α13*Abs[A3[x]]^2 + α14*Abs[A4[x]]^2), 
          A2'[x] == I*κ2*Conjugate[A1[x]]*A3[x]*A4[x]*E^(-I*Δkl*x) + I*A2[x]*(α21*Abs[A1[x]]^2 + α22*Abs[A2[x]]^2 + α23*Abs[A3[x]]^2 + α24*Abs[A4[x]]^2), 
          A3'[x] == I*κ3*A1[x]*A2[x]*Conjugate[A4[x]]*E^(I*Δkl*x) + I*A3[x]*(α31*Abs[A1[x]]^2 + α32*Abs[A2[x]]^2 + α33*Abs[A3[x]]^2 + α34*Abs[A4[x]]^2), 
          A4'[x] == I*κ4*A1[x]*A2[x]*Conjugate[A3[x]]*E^(I*Δkl*x) + I*A4[x]*(α41*Abs[A1[x]]^2 + α42*Abs[A2[x]]^2 + α43*Abs[A3[x]]^2 + α44*Abs[A4[x]]^2), 
          A1[0] == (I0*25)/ω1, A2[0] == (I0*25)/ω2, A3[0] == 0, A4[0] == j};

DEsols = ParametricNDSolve[system, {A1[x], A2[x], A3[x], A4[x]}, {x, 0, 2000}, {j}]
Plot[Evaluate@Table[Abs[(A4[j][x]) /. DEsols]^2, {j, 0, 10}], {x, 0, 2000}, PlotStyle -> {Orange}, PlotLegends -> {"A4"}, PlotRange -> All, AxesOrigin -> {0, 0}]

However, it is not plotting and I'm not sure what I've done wrong. Furthermore, I intend to plot A4[x] as a function of A4[0] for a fixed x (x=2000). How should I go about fixing this? Thank you.

Answer 1

After modifying ParametricNDSolve to ParametricNDSolveValue (without argument brackets ...[x])

DEsols = ParametricNDSolveValue[system, {A1 , A2 , A3 , A4 }, {x, 0, 2000}, {j}]

you can access A4

Plot[sols[0][[4]][x], {x, 0, 2000}]

Unfortunately your system seems to be solvable only for j==0

Answer 2

We can use NDSolve[] to solve this problem as follows:

system[j_] := {A1'[x] == 
    I*\[Kappa]1*Conjugate[A2[x]]*A3[x]*A4[x]*
      E^(-I*\[CapitalDelta]kl*x) + 
     I*A1[x]*(\[Alpha]11*Abs[A1[x]]^2 + \[Alpha]12*
         Abs[A2[x]]^2 + \[Alpha]13*Abs[A3[x]]^2 + \[Alpha]14*
         Abs[A4[x]]^2), 
   A2'[x] == 
    I*\[Kappa]2*Conjugate[A1[x]]*A3[x]*A4[x]*
      E^(-I*\[CapitalDelta]kl*x) + 
     I*A2[x]*(\[Alpha]21*Abs[A1[x]]^2 + \[Alpha]22*
         Abs[A2[x]]^2 + \[Alpha]23*Abs[A3[x]]^2 + \[Alpha]24*
         Abs[A4[x]]^2), 
   A3'[x] == 
    I*\[Kappa]3*A1[x]*A2[x]*Conjugate[A4[x]]*
      E^(I*\[CapitalDelta]kl*x) + 
     I*A3[x]*(\[Alpha]31*Abs[A1[x]]^2 + \[Alpha]32*
         Abs[A2[x]]^2 + \[Alpha]33*Abs[A3[x]]^2 + \[Alpha]34*
         Abs[A4[x]]^2), 
   A4'[x] == 
    I*\[Kappa]4*A1[x]*A2[x]*Conjugate[A3[x]]*
      E^(I*\[CapitalDelta]kl*x) + 
     I*A4[x]*(\[Alpha]41*Abs[A1[x]]^2 + \[Alpha]42*
         Abs[A2[x]]^2 + \[Alpha]43*Abs[A3[x]]^2 + \[Alpha]44*
         Abs[A4[x]]^2), A1[0] == (I0*25)/\[Omega]1, 
   A2[0] == (I0*25)/\[Omega]2, A3[0] == 0, A4[0] == j}; 
sol[j_] := NDSolve[system[j], {A1, A2, A3, A4}, {x, 0, 2000}]

Now we can plot solution for small range of j

 {LogLogPlot[
  Evaluate@Table[
    Abs[(A1[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A1"], 
 LogLogPlot[
  Evaluate@Table[
    Abs[(A2[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A2"], 
 LogLogPlot[
  Evaluate@Table[
    Abs[(A3[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A3"], 
 LogLogPlot[
  Evaluate@Table[
    Abs[(A4[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A4"]}