# A fraught with incorrect results ODE

I mean the following ODE $$y''(x)+y'(x)=\exp (-2 x) y(x)^3.$$

Trying to solve it in version 13.1 on Windows 10 by

DSolve[y''[x] + y'[x] == Exp[-2 x]*y[x]^3, y[x], x]


, I obtain a huge incorrect result

{{y[x] -> 1/4 (-(E^( 1/2 (C - Inactive[Integrate][(E^(-2 K) (y[K]^5 + 2 C y[K]^5 + 2 y[K]^5 Inactive[Integrate][( E^(-2 K) (y[K] - Derivative[y][K]) (-y[K]^4 + E^(2 K) y[K] Derivative[y][K] + E^(2 K) Derivative[y][K]^2))/( 2 y[K]^3), {K, 1, K}] + 8 E^(2 K) C^2 y[K]^3 Sqrt[ 1 + 4 C + ...

Pay your attention to Inactive[ Integrate][(E^(-2 K) (y[K] - Derivative[y][K]) (-y[K]^4 + E^(2 K) y[K] Derivative[y][K] + E^(2 K) Derivative[y][K]^2))/(2 y[K]^3), {K, 1, K}]

in the above, where the function y[x] is expressed through itself and its derivative y'[x].

Next, the command

DSolve[{y''[x] + y'[x] == Exp[-2 x]*y[x]^3, y == 1, y' == -1},  y[x], x]


is running without any response for hours. Likely an infinite loop is created since the resourсes of my comp are not exhausted.

The change of the independent variable x by

DSolveChangeVariables[ Inactive[DSolve][{y''[x] + y'[x] == Exp[-2 x]*y[x]^3, y == 1,
y' == -1}, y[x], x], u, t, t == Exp[-x]]


Inactive[DSolve][{t u[t]^3 == t u''[t], u == 1, DSolve'DSolveChangeVariablesDump'd$18576[u] == -1}, u[t], t] produces at least two bugs: u==1 instead of u==1 and DSolve'DSolveChangeVariablesDump'd$18576[u] == -1.

The questions arise: how to correctly solve this ODE? are there workarounds?

As it has been noticed in the original question one can transform the independent variable $$x \to t=\exp(-x)$$. Then our differential equation transforms to $$\;u''(t)=u(t)^3$$ and the initial conditions transforms to $$\;u(1)=1\;$$ and $$\;u'(1)=1$$, where $$u(t)=\tilde{y}(\exp(-x))=y(x)$$. Physicists usually oversimplify notation as e.g. $$\;y(x)= y(\exp(-x))=y(t)$$, however we should carefully distinguish all the functions $$u, \tilde{y}, y$$. Derivation of the transformed equation is a simple excersise for users using version 13.0 and earlier.

The transformed equation can be integrated once multiplying it by $$u'(t)$$: $$0=u''(t) - u(t)^3=u' u'' -u^3 u'=(\frac{1}{2} {u'}^2-\frac{1}{4}u^4+c)'=0$$ where $$c$$ is a constant. Solution to this equation can be reformulated as integration of an elliptic integral see e.g. Solving 0=−λ ϕ(t)^3+μ^2 ϕ(t)+ϕ′′(t) and usually Mathematica solves such equations this way obtaining involved solutions in terms of inverse functions. We wouldn't like solving it this way and we have to play with further transforming it to another form which the system can recognize. Usually I prefer transforming to the canonical Weierstrass form (see e.g. 1, 2, 3, 4, ...) then we encounter problems with solving differential equations with given initial (or boundary) conditions see e.g. 5, nevertheless there are still appropriate tools one can harness e.g. 6. Such a procedure allows us to get exact solutions with a certain number of intermediate steps, reflecting various problematic issues when solving differential equations in terms of elliptic functions.

Nonetheless it appears that we can solve our initial value problem directly prescribing initial conditions

us[t_] = FullSimplify @ DSolveValue[{u''[t] - u[t]^3 == 0, u == 1, u' == 1},
u[t], t]

-(-1)^(1/4) JacobiSN[(-(1/2) + I/2) (-1 + t)
+ InverseJacobiSN[(-1)^(3/4), -1], -1]


This is a real function over a real domain:

Plot[ Flatten[{ReIm[-(-1)^(1/4) JacobiSN[(-(1/2) + I/2) (-1 + t) +
InverseJacobiSN[(-1)^(3/4), -1], -1]], t}],
{t, 0, 1.6}, Evaluated -> True, AspectRatio -> Automatic,
Epilog -> {Red, PointSize[0.02], Point[{1, 1}]}] being an elliptic function (doubly periodic meromorphic function in the complex plane see e.g. 6, 7) however the solution to the oringinal equation is composed with an exponential function and so it isn't an elliptic function anymore:

ys[x_] = us[t] /. t -> Exp[-x]


it asymptotically goes at infinity to

us // N // Chop

0.219982

Plot[ Flatten[{ ReIm[ys[x]], Exp[-x], us // Chop}], {x, -1/2, 3},
Evaluated -> True, AspectRatio -> Automatic] Edit

Another question arised in the comments, namely how we can find appropriate parameters in general solution found with DSolve to be compatibile with exact solution to the initial value problem found above i.e. us[t]? In general we need not play with the full symbolic power of the system but we would rather take a shortcut approach solving numerically appropriate equations finding parameters. This is an analogous problem to e.g. How to remove irrelevant terms (such as log 's) in the solution of differential equation?, however here instead of elementary transcendental functions as Log we have to deal with higher transcendental functions as JacobiSN what makes the problem more difficult.

Let's define the first general solution to $$\;u''(t)=u(t)^3$$:

u1[t_, c1_, c2_] =
u[t] /. First @ FullSimplify @ DSolve[u''[t] - u[t]^3 == 0, u[t], t]/.
{C -> c1, C -> c2} //Quiet

-((I 2^(1/4) JacobiSN[-(((1 - I) Sqrt[Sqrt[c1] (c2 + t)^2])/2^(
3/4)), -1])/Sqrt[(I/Sqrt[c1])])


Since we have two parameters we should compare u1 and us for two different arguments providing non-degenerate system of equations choosing appropriate starting points in FindRoot, e.g.

{c1c, c2c} = {c1, c2} /. FindRoot[{u1 == us, u1[1/2] == us[1/2]},
{{c1, 1}, {c2, 0}}] // Chop

{0.5, 0.311029}


now we can see that the both solutions coincide numerically:

Plot[{us[t], Re @ u1[t, c1c, c2c]}, {t, 0, 3/2},
PlotStyle -> {{Thick}, {Dashed, Green}}, PlotLegends -> "Expressions"] A symbolic approach would involve searching for appropriate formulas in e.g. Entity["MathematicalFunction", "JacobiSN"]["Dataset"] and show that the both functions are equal under specific conditions but that is a different story. Here we can simply identify c1 as 1/2 and comparing u1[t, 1/2, c2] with us[t]} we find also an exact formula for c2:

c2 /. First @ Solve[ (-(1/2) + I/2) (-1) + InverseJacobiSN[(-1)^(3/4), -1]
== (-(1/2) + I/2) c2, c2]

 -1 - (1 + I) InverseJacobiSN[(-1)^(3/4), -1]

• How to derive it from {{{u[E^-x] -> -2^(1/4) Sqrt[I/Sqrt[C]] Sqrt[C] JacobiSN[((-1)^(3/4) Sqrt[ Sqrt E^(-2 x) Sqrt[C] + 2 Sqrt E^-x Sqrt[C] C + Sqrt Sqrt[C] C^2])/Sqrt, -1]}, {u[E^-x] -> 2^(1/4) Sqrt[I/Sqrt[C]] Sqrt[C] JacobiSN[((-1)^(3/4) Sqrt[ Sqrt E^(-2 x) Sqrt[C] + 2 Sqrt E^-x Sqrt[C] C + Sqrt Sqrt[C] C^2])/Sqrt, -1]}}? Oct 23, 2022 at 7:27
• Or {{u[t] -> -2^(1/4) Sqrt[I/Sqrt[C]] Sqrt[C] JacobiSN[((-1)^(3/4) Sqrt[ Sqrt t^2 Sqrt[C] + 2 Sqrt t Sqrt[C] C + Sqrt Sqrt[C] C^2])/Sqrt, -1]}, {u[t] -> 2^(1/4) Sqrt[I/Sqrt[C]] Sqrt[C] JacobiSN[((-1)^(3/4) Sqrt[ Sqrt t^2 Sqrt[C] + 2 Sqrt t Sqrt[C] C + Sqrt Sqrt[C] C^2])/Sqrt, -1]}}? Oct 23, 2022 at 7:29
• @user64494 This is a general solution of the underlying ODE without imposed initial conditions. If you have the above solution ys[x] you could find parameters C, C numerically with e.g. FindRoot. Symbolically I believe there is no an automatic approach. Ideally one could apply Reduce with appropriate conditions, however I'm almost sure that currently it cannot succeed. You should rather look for appropriate formulas among Entity["MathematicalFunction", "JacobiSN"]["Dataset"] e.g. in "ArgumentSimplifications". Oct 23, 2022 at 13:55
• Artes (@ does not work.): You wrote "If you have the above solution ys[x] you could find parameters C1], C numerically with e.g. FindRoot". I unsuccessfully tried it with Mathematica and Maple. Oct 23, 2022 at 14:14
• @user64494 Analogous issue came up e.g. here where FindRoot suceeded while Reduce did not. At the moment I have no time to dmemonstrate it but I guess I could do it perhaps with a slightly more involved approach. Oct 23, 2022 at 14:40

Mathematica V 13.1 can not solve it. (well, it solved it, but the solution it gave, which you showed, is not really useful).

May be next version will. I did not see a workaround in Mathematica.

Here is the solution. I converted it from Maple which can solve it.

It was little tricky to do the conversion since Maple uses little different definition of the JacobiSN function from Mathematica and took me a little while to find the difference so it works OK in Mathematica and satisfies the ode. The difference turned out minor but not documented. Any way, here is the solution in Mathematica syntax, which verifies OK

ClearAll[y, x];
ode = y''[x] + y'[x] - Exp[-2*x]*y[x]^3 == 0
sol = y -> Function[{x}, C*JacobiSN[(-1/2*Sqrt[-2*Exp[-2*x]] + C)*C, -1]]
ode /. sol // FullSimplify

(* True *)


Solution in Latex

$$c_2 \operatorname{JacobiSN}\left(\left(c_1-\frac{\sqrt{-e^{-2 x}}}{\sqrt{2}}\right) c_2,-1\right)$$

If you want the actual method/steps used to solve the ode, then I think the math forum will be the right place in this case. Maple does not give a hint on how it solved it, other than saying it tried JacobiSN and it "worked".

• You wrote "well, it solved it, but the solution it gave, which you showed, is not really useful)". Let us call things by their proper names: that solution is extremely wrong and causes serious doubts about reliability of some Mathematica results. Oct 19, 2022 at 3:28
• Thank you for your comment. I was aware about the Maple's solution before having submitted the question. If I correctly understand it, the result of infolevel[all] := 4; dsolve(diff(y(x), x, x)+diff(y(x), x) = exp(-2*x)*y(x)^3, y(x))  suggests Maple uses differential algebra to this end. It should be noticed I can derive any non-zero particular solution from the general one neither in Mathematica nor in Maple. Oct 19, 2022 at 3:31
• For example, I mean the particular solution for {y''[x] + y'[x] == Exp[-2 x]*y[x]^3, y == 1, y' == -1}. Oct 19, 2022 at 4:25

Mathematica 13.1 is able to solve that ODE after the change of the independent variable x.

DSolveChangeVariables[
Inactive[DSolve][y''[x] + y'[x] == Exp[-2x] y[x]^3, y[x], x], u, t, t == Exp[-x]]

Inactive[DSolve][t u[t]^3 == t u''[t], u[t], t]

Activate[Inactive[DSolve][t u[t]^3 == t u''[t], u[t], t],
Unevaluated[DSolve]] /. t -> Exp[-x]

{{{u[E^-x] -> -2^(1/4) Sqrt[I/Sqrt[C]] Sqrt[C]
JacobiSN[((-1)^(3/4) Sqrt[Sqrt E^(-2 x) Sqrt[C] +
2 Sqrt E^-x Sqrt[C] C + Sqrt Sqrt[C] C^2])/Sqrt, -1]},
{u[E^-x] -> 2^(1/4) Sqrt[I/Sqrt[C]] Sqrt[C]
JacobiSN[((-1)^(3/4) Sqrt[Sqrt E^(-2 x) Sqrt[C] +
2 Sqrt E^-x Sqrt[C] C + Sqrt Sqrt[C] C^2])/Sqrt, -1]}}


Unfortunately, I cannot derive the particular solution for {y''[x] + y'[x] == Exp[-2 x] y[x]^3, y == 1, y' == -1} from it.

While for this example, the output is exceedingly complicated and unhelpful sometimes the best one can hope for is an implicit solution. As a simple example (from the documentation), evaluate the following example to see the implicit solution to an Abel equation:

DSolve[y'[x] == y[x]^3 - ((x + 1) y[x]^2)/x, y[x], x] // Quiet


Back to your ODE: Twenty years ago I co-authored a paper on The Jacobi elliptic-function method for finding periodic wave solutions to nonlinear evolution equations. The code for that paper—my most cited work—was written in Mathematica and appears to run fine in the latest version. And the method used in that paper should quickly find solutions to ODEs such as this one (possibly still requiring a natural change of variables).

Aside: I'd be interested to know how your ODE arose.

• I saw it at a math forum in Russian where the solution by hand was discussed for a special case. Nov 26, 2022 at 11:51