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

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.

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

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


• Thanks for the response. However when I try to plot A4[2000] as a function of j, a bunch of warnings and errors showed up. In particular, I did LogLinearPlot[Abs[(A4[2000]) /. sol[j]]^2, {j, 0, 50}, PlotRange -> All] and it refuses to plot. Commented Jul 9, 2020 at 16:35
• @kowalski Just check that my plots made for {j, .5 10^-10, 2 10^-10, .5 10^-10}. And for this small parameter we have A4 of order 10^29. With increasing j over 3*10^-6 there is no stable solution. Commented Jul 9, 2020 at 16:59
• I can see from the fourth panel A4 that it is indeed in the order of 10^29. However, if I specify my methods, namely sol[j_] := NDSolve[system[j], {A1, A2, A3, A4}, {x, 0, 2000}, Method -> "Extrapolation", StartingStepSize -> 0.01] over j up to 3*10^-6 I do see a plot but I am unsure if it's the right plot since there are many warnings Commented Jul 9, 2020 at 17:22
• @kowalski What actually do you try to describe by this system? Commented Jul 9, 2020 at 18:18
• A1, A2, A3, A4 are amplitudes of a wave in an optic fiber amplifier. A1 and A2 represents the pump amplitudes, A3 represents the idler amplitude and A4 represents the signal amplitude. I am trying to plot gain which is Log10[Abs[A4[x]/A4[0]]^2] as a function of A4[0] which is j to demonstrate saturation of gain as a function of input signal power (Abs[A4[0]]^2) Commented Jul 9, 2020 at 18:23