Tag: 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.