# NDSolve stiff system of ODEs

I am having a hard time trying to solve the following system of ODEs

X'' = e^{2X} + e^{2Y}
Y'' = e^{2Y} - 3 e^{2X}


with initial conditions that mimic a linear behavior of the solutions at $$-\infty$$. I am using for instance something like

X(u0) = lambdaX * u0
Y(u0) = lambdaY * u0
X'(u0) = lambdaX
Y'(u0) = lambdaY


with $$u0 = -10$$. The actual values of the parameters $$\lambda_X$$ and $$\lambda_Y$$ do not matter much, as long as the solutions are linear at $$-\infty$$ (ie. both the exponentials have to go to $$0$$ at $$-\infty$$, so that $$X'' \rightarrow 0$$ and $$Y'' \rightarrow 0$$). I tried for instance with $$\lambda_Y = 1.5$$ and $$\lambda_X = 0.5$$.

The basic code that I'm using is something like:

lambdaY = 1.5;
lambdaX = 0.5;
u0 = -10;
sol = NDSolve[{X''[u] == Exp[2 Y[u]] + Exp[2 X[u]], Y''[u] == Exp[2 Y[u]] - 3 Exp[2 X[u]], Y[u0] == u0*lambdaY, X[u0] == u0*lambdaX, Y'[u0] == lambdaY, X'[u0] == lambdaX}, {Y, X}, {u, u0, 10}, Method -> "StiffnessSwitching", AccuracyGoal -> \[Infinity], MaxSteps -> 5000000, InterpolationOrder -> All]


The problem is the usual step size is effectively zero; singularity or stiff system suspected.: the solutions are diverging suddenly at some point and I could not understand what is happening there. I have tried with the standard ways to fix the issue, ie. playing with WorkingPrecision, AccuracyGoal and PrecisionGoal and changing the Method to StiffnessSwitching but none of these worked. I also tried to redefine the variables like

R = e^{2X}
S = e^{2Y}


but I did not manage to avoid the problem.

I am looking for a solution that is asymptotically linear at $$-\infty$$ for both functions $$X$$ and $$Y$$. I am sure that the solution is a simple one and I am missing something. Any suggestion on how to solve this problem?

Thanks a lot.

• Be careful with the use of the _ character that is shortcut for Blank, which is a is a pattern object that can stand for any Wolfram Language expression, and consequently, not suitable for variable names. Commented May 15, 2023 at 12:07
• Thanks for the comment! I don't think this is the problem though, in the actual Mathematica code I'm not using underscores.
– FM87
Commented May 15, 2023 at 12:48
• You probably have a singularity, not stiffness. Can't confirm without working code, though. Commented May 15, 2023 at 14:02
• Thanks for the help. I added a working code snipped in the original post.
– FM87
Commented May 15, 2023 at 14:19

If we first consider $$X''(u)=e^{2 X(u)} \,,$$ which has a solution $$\log \left(\frac{1}{2} c_1 \left| \sec \left(\frac{1}{2} c_1 (u-{c_2})\right)\right| \right) \,,$$ we see it approaches infinity in finite time.

In the OP's system, $$X''(u)=e^{2 X(u)}+e^{2 Y(u)} > e^{2 X(u)} \,.$$ We should expect that once $$X(u)$$ starts increasing, it shortly goes to infinity, too.

This can be seen in the plots of $$X$$, $$Y$$ and their derivatives:

GraphicsGrid@Map[
ListLinePlot[# /. First@sol, PlotLabel -> #[u], Frame -> True] &,
{{X, Y}, {X', Y'}, {X'', Y''}}, {2}]


The orders of growth of $$X$$, $$X'$$, $$X''$$ near the singularity are consistent with a logarithmic pole, too, similar to the solution of the ODE example that began this answer.

Here's a way to see the orders of growth at the singularities $$u=u_2 \approx -0.0018$$:

With[{u2 = Last@First[X@"Domain" /. First@sol]},
Function[{xy, scale},
Plot[scale #[u] /. First@sol, {u, -0.05, u2}, PlotLabel -> #[u],
Frame -> True] & /@ xy],
{{{X, Y}, {X', Y'}, {X'', Y''}},
{1/Log[u2 - u], u2 - u, (u2 - u)^2}}]
]


$$X$$, $$X'$$, $$X''$$ are asymptotic to $$-\log(u_2-u)$$, $$1/(u_2-u)$$, $$1/(u_2-u)^2$$ respectively; $$Y$$, $$Y'$$, $$Y''$$ are asymptotic to $$3\log (u_2-u)$$, $$-3/(u_2-u)$$, $$-3/(u_2-u)^2$$ respectively.

• Thanks, very good point! I had not realized that actually the simpler equation X'' = Exp[2X] diverges at finite u. That explains a divergence also in the equation I'm looking at.
– FM87
Commented May 15, 2023 at 22:45
• @FM87 You're welcome. Commented May 15, 2023 at 22:47