# Two first order non-linear differential equation: Using DSolve to obtain the solutions analytically

I have two simple first order (coupled) differential equations which I denote as system together with its initial conditions

system = {A1'[x] == (j1/2)*A1[x]^2*Conjugate[A1[x]] + j1*A1[x]*A2[x]*Conjugate[A2[x]],
A2'[x] == (j2/2)*A2[x]^2*Conjugate[A2[x]] + j2*A1[x]*A2[x]*Conjugate[A1[x]],
A1 == √[P10], A2 == √[P20];


where j1, j2, P10, P20 are all constants and independent of x. I then do

DSolve[system,{A1,A2},x]


And I'm returned with: There are fewer dependent variables than equations, so the system is overdetermined. But I have exactly 2 equations with 2 unknowns, so I'm not sure why Mathematica is returning this error. What am I doing wrong here?

Thank you.

• Is \Sqrt a typo or you really write this in your code? – xzczd Sep 7 at 2:06
• Which variables and constants, if any, are complex? Do you have reason to suppose that analytical solutions actually exist? DSolve is unable to find a solution, although it does not generate an error message for me. – bbgodfrey Sep 7 at 3:56
• Try quitting the kernel and then retry your code. – Michael E2 Sep 7 at 5:38
• \Sqrt is a typo, but I didn't type that wrong in my code @xzczd – kowalski Sep 7 at 6:29
• I mentioned j1, j2, P10 and P20 are constants @bbgodfrey An analytical solution should exist given that it represents the pump amplitude of a physical system. Even if it doesn't it should return something that says unable to find solution. Not what im getting. Also quitting the kernel didn't change anything – kowalski Sep 7 at 6:32

Amazingly, this system can be solved symbolically, at least for j2 = j1. How useful the answer may be is a different issue. If all symbols are real, then system (without initial conditions) can be rewritten as

system = {p1'[x] == j1 p1[x]^2 + 2 j1 p1[x] p2[x],
p2'[x] == 2 j2 p1[x] p2[x] + j2 p2[x]^2}


where p1[x] = A1[x]^2 and p2[x] = A2[x]^2. (If the symbols are complex, then p1[x] = A1[x]*Conjugate[A1[x]], p2[x] = A2[x]*Conjugate[A2[x]], and j1 and j2 replaced by their real parts.) The obvious next step, DSolve[system /. j2 -> j1, {p1, p2}, x] returns unevaluated with unexpected error messages. So, instead, eliminate p2, which is straightforward.

D[Simplify[#/p1[x]], x] & /@ First[system];
% /. Rule @@ Last[system];
Solve[First[system], p2[x]] // Flatten;
eq = Simplify[%% /. %]
(* 3 j1 j2 p1[x]^3 + 2 p1''[x] ==
2 (j1 + j2) p1[x] p1'[x] + ((2 j1 + j2) p1'[x]^2)/(j1 p1[x]) *)


which can be solved.

s = DSolve[eq /. j2 -> j1, b1[x], x][[1, 1]]
(* p1[x] -> InverseFunction[Inactive[Integrate][4/(3 (4 j1 K^2 +
(E^((3 C)/8) j1^2 K^3)/(4 E^((3 C)/8) j1^3 K^5 +
Sqrt[-E^(((9 C)/8)) j1^6 K^9 + 16 E^((3 C)/4) j1^6 K^10])^(1/3) +
(4 E^((3 C)/8) j1^3 K^5 + Sqrt[-E^(((9 C)/8)) j1^6 K^9 +
16 E^((3 C)/4) j1^6 K^10])^(1/3))), {K, 1, #1}] &][x + C] *)


Perhaps, more convenient is

s[[2, 1]] - C -> s[[2, 0, 1]][s[]] - C;
% /. {K, 1, p1[x]} -> {K, p10, p1[x]} /. C -> 0
(* x -> Inactive[Integrate][4/(3 (4 j1 K^2 +
(E^((3 C)/8) j1^2 K^3)/(4 E^((3 C)/8) j1^3 K^5 +
Sqrt[-E^(((9 C)/8)) j1^6 K^9 + 16 E^((3 C)/4) j1^6 K^10])^(1/3) +
(4 E^((3 C)/8) j1^3 K^5 + Sqrt[-E^(((9 C)/8)) j1^6 K^9 +
16 E^((3 C)/4) j1^6 K^10])^(1/3))), {K, p10, p1[x]}] *)


p2 then can be determined by inserting the p1solution into system[], and the initial condition p20 applied to obtain C. If j2!= j1, DSolve returns unevaluated.

C can be derived more simply by obtaining a first integral of eq by the standard transformation valid for autonomous second order ODEs. Doing so provides some useful insight as well.

eqv = eq /. {p1''[x] -> v'[z] v[z], p1'[x] -> v[z], p1[x] -> z}
(* 3 j1 j2 z^3 + 2 v[z] v'[z] == (j1 + j2) z v[z] + ((2 j1 + j2) v[z]^2)/(j1 z) *)

eqc1 = (DSolve[eqv, v[z], z] // First // Simplify) /. {v[z] -> p1'[x], z -> p1[x]}
(* (1/(3 j2)) 2 (j1 + j2)^2 (Log[p1[x]]/j1 +
(Log[(2 (j1 + j2) (j1 p1[x]^2 - p1'[x]))/(j1 p1[x]^2)] - (3 j2
Log[-((2 (j1 + j2) (3 j1 j2 p1[x]^2 + (-2 j1 + j2) p1'[x]))/(3 j1 j2 p1[x]^2))])
/(2 j1 - j2))/(j1 - 2 j2)) == C *)


from which C is obtained by setting x == 0 and applying initial conditions. To obtain the complete solution for p1[x] requires solving eqc1 for p1'[x], which can be done for j2 = j1 and, perhaps, a few other cases. (This explains why DSolve[eq, b1[x], x] is unable to obtain solutions for arbitrary choices of j2/j1.) So, consider j2 = j1.

Simplify[eqc1 /. {p1'[x] -> j1 p10^2 + 2 j1 p10  p20, p1[x] -> p10}]
Simplify[% /. j2 -> j1]
Exp[3/8 #] & /@ %
(* -((64 j1^2 (-p10 + p20)^3)/(27 p10 p20)) == E^((3 C)/8) *)


This procedure also provides a simple evaluation of p2[x]. Interchange p1[x] and p2[x], and p10 and p20 in the code block just before this addendum and also in the expression for C just above.