# Why does Mathematica give no output?

I am trying to solve a system of coupled differential equations. I am using the following code:

DSolve[
{D[A[z], z] + B[z]*(Λ - Sqrt[2]*b) - (a/Sqrt[2])*(A[z]^2 - B[z]^2)*B[z] == 0,
D[B[z], z] - A[z]*(Λ + Sqrt[2]*b) - (a/Sqrt[2])*(A[z]^2 - B[z]^2)*A[z] == 0},
{A[z], B[z]}, z]


Upon executing, Mathematica gives no output at all. It does not even give an error message. What's wrong?

• Mathematica 11.2 returns unevaluated, probably because it cannot solve this problem. Many nonlinear ODEs have no solution. – bbgodfrey Dec 25 '17 at 16:34
• putting assumptions on your parameters (real, positive) might help, but not likely – george2079 Dec 25 '17 at 16:41

Because this system of ODEs is autonomous, z can be eliminated, reducing the system to a single first order ODE.

eq1 = D[A[z], z] == -B[z]*(Λ - Sqrt[2]*b) + (a/Sqrt[2])*(A[z]^2 - B[z]^2)*B[z];
eq2 = D[B[z], z] == +A[z]*(Λ + Sqrt[2]*b) + (a/Sqrt[2])*(A[z]^2 - B[z]^2)*A[z];


Specifically, replace A[z] by x and B[z] by y[x]:

eq = eq2 /. B'[z] -> y'[x] A'[z] /. A'[z] -> eq1[[2]] /. {B[z] -> y[x], A[z] -> x}

(* (-(-Sqrt[2] b + Λ) y[x] + (a y[x] (x^2 - y[x]^2))/Sqrt[2]) y'[x] ==
x (Sqrt[2] b + Λ) + (a x (x^2 - y[x]^2))/Sqrt[2] *)


from which DSolve obtains a parametric solution.

sol = Values@DSolve[eq, y[x], x] // Flatten

{-Sqrt[(2 b)/a + x^2 - (Sqrt[2] Λ)/a - (Sqrt[2] Sqrt[2 b^2 - 2 Sqrt[2] b Λ -
2 Sqrt[2] a x^2 Λ + Λ^2 + Sqrt[2] a C[1]])/a],
Sqrt[(2 b)/a + x^2 - (Sqrt[2] Λ)/a - (Sqrt[2] Sqrt[2 b^2 - 2 Sqrt[2] b Λ -
2 Sqrt[2] a x^2 Λ + Λ^2 + Sqrt[2] a C[1]])/a],
-Sqrt[(2 b)/a + x^2 - (Sqrt[2] Λ)/a + (Sqrt[2] Sqrt[2 b^2 - 2 Sqrt[2] b Λ -
2 Sqrt[2] a x^2 Λ + Λ^2 + Sqrt[2] a C[1]])/a],
Sqrt[(2 b)/a + x^2 - (Sqrt[2] Λ)/a + (Sqrt[2] Sqrt[2 b^2 - 2 Sqrt[2] b Λ -
2 Sqrt[2] a x^2 Λ + Λ^2 + Sqrt[2] a C[1]])/a]}


Here is a sample plot

Plot[sol /. {Λ -> 1, a -> 1, b -> 1} /. C[1] -> 2, {x, -1.1, 1.1}, AxesLabel -> {A, B},
ImageSize -> Large, LabelStyle -> Directive[Bold, Black, Medium], AspectRatio -> Full]


This result can be validated by direct numerical solution of the original equations.

s = ParametricNDSolveValue[{eq1, eq2, A[0] == A0, B[0] == B0}, {A, B}, {z, 0, 10},
{A0, B0, a, b, Λ}];

params = Sequence[1, 1, 1, 1, 1];
Plot[Evaluate@Through[s[params][z]], {z, 0, 10}, AxesLabel -> {z, "A, B"},
ImageSize -> Large, LabelStyle -> Directive[Bold, Black, Medium]]
ParametricPlot[Evaluate@Through[s[params][z]], {z, 0, 10}, AxesLabel -> {A, B},
ImageSize -> Large, LabelStyle -> Directive[Bold, Black, Medium], AspectRatio -> Full]


and a parametric plot identical to the first plot in this answer.

Given that sol provides B as a function of A, in principle A can be determined as a function of z from

Simplify[eq1 /. B[z] -> sol[[j]] /. x -> A[z]]


where j represents each of the four components of sol. For instance, with j == 1,

(* Sqrt[2 b^2 - 2 Sqrt[2] b Λ + Λ^2 - 2 Sqrt[2] a Λ A[z]^2 + Sqrt[2] a C[1]]
Sqrt[(2 b + a A[z]^2 - Sqrt[2] (Λ + Sqrt[2 b^2 - 2 Sqrt[2] b Λ + Λ^2 -
2 Sqrt[2] a Λ A[z]^2 + Sqrt[2] a C[1]]))/a] + A'[z] == 0 *)


which can be solved with DSolveValue or Integrate. Unfortunately, neither gives a closed-form solution for A[z]. In summary, this answer derives a parametric relation between A and B, but does not provide closed-form solutions for A or B as functions of z. Whether the parametric solution is sufficient depends on the OP's goal.