# Non-linear system of differential equation with stiff potential

I am trying a system of two differential equations that looks pretty simple, but with a potential varying a lot on a short period centered on 0. The definition of the potential followed by the attempt at solving the differential equation is

a0[y_] := 10^(-22) y^2/k^2 + Sqrt[1 + y^2/k^2];
a2[y_] := D[D[a0[y], y], y];
f[y_] := 10^(-18)* a2[y]/a0[y]^3;
g[y_] := 1 + f[y];
pot[y_] := D[D [g[y], y], y]/g[y]; (* Potential *)

yi = -300;
yf = 100;
k = 10^(-20);
sol = NDSolve[{x1'[y] == (pot[y] - 1) - x1[y]^2 + x2[y]^2,
x2'[y] == -2*x1[y]*x2[y], x1[yi] == 0, x2[yi] == -1} , {x1,x2}, {y, yi, yf}, WorkingPrecision -> 30, MaxSteps -> Infinity][[1]]
Plot[Evaluate[x1[y] /. sol], {y, yi, yf}, PlotRange -> All, PlotPoints -> 500, AxesLabel -> {y, "x1(y)"}]
Plot[Evaluate[x2[y] /. sol], {y, yi, yf}, PlotRange -> All, PlotPoints -> 500, AxesLabel -> {y, "x2(y)"}]


And it appears that the solutions $$x1$$ and $$x2$$ are just equal to the initial conditions. I was expecting something growing after 0, since the potential there is about $$10^{40}$$, but I can't even see a little deviation from the initial values. Any comment, any clue on what I am doing wrong is greatly appreciated!

• I think you need to transform your differential equation (substitutions) in order to get closer to a numerically solvable model. Jun 3, 2019 at 11:30
• If you change a2[y]:=... to a2[y_] := Derivative[2][a0][y]; and pot[y]:=... to pot[y_] := Derivative[2][g][y]/g[y] the potential can be evaluated and evaluates to Ò[10^-68]. That might explain the behavior. Jun 3, 2019 at 11:36
• What I mean is that your pot[y] is so singular at y=0 that you could simplify your differential equation dramatically by approximating pot[y_] = q*DiracDelta[y] (for some amplitude q) and trying to find an analytic solution of the differential equation. This analytic solution could then be used as a base to find the exact solution by appropriate transformations. Jun 3, 2019 at 13:13
• @Roman I tried what you suggested and, unsurprisingly, NDSolve encountered a step size effectively zero (for y around pi). I then tried to solve the system using DSolve with and without initial conditions, but the output MMA gave me was the input I gave. I am puzzled here... Jun 4, 2019 at 14:49
• @Free_ion for the Dirac $\delta$-suggestion you could solve the equations analytically, it may be easier than NDSolve. Are you familiar with the techniques for doing this in the presence of $\delta$-potentials? Essentially you solve the differential equation for $pot=0$ and then connect two different regions ($y<0$ and $y>0$) by a kink given by the $\delta$-function amplitude. Jun 4, 2019 at 14:58

Just to illustrate my comment:

Clear[a2, pot]
a2[y_] := Derivative[2][a0][y];
pot[y_] := Derivative[2][g][y]/g[y];(*Potential*)

Plot[pot[y], {y, yi, yf}, PlotRange -> All]


and

Plot[pot[y], {y, yi, yf}, PlotRange -> All,PlotPoints -> {Automatic, {0}}]


• Well, the potential itself is well-defined, and is actually similar to this one link, with different ranges in x and y. But changing the definition of the derivative does not give a different behaviour for x1 and x2, unfortunately. Jun 4, 2019 at 12:34
• The change from D[..] to Derivative[..] changes a lot . Look at pot[y] in both versions. Only the second one can be evaluated/plotted! Jun 4, 2019 at 12:43
• But is this really an effect of using Derivative[..] instead of D[..]? It seems to me the only difference is in the number of points used to make the plot. Jun 4, 2019 at 12:55
• Try to evaluate  pot[0] and you'll get, expecting a number, a strange result \!$$\*SubscriptBox[\(\[PartialD]$$, $$0$$]$$\*SubscriptBox[\(\[PartialD]$$, $$0$$]$$(1 + \*FractionBox[\( \*SubscriptBox[\(\[PartialD]$$, $$0$$]$$\*SubscriptBox[\(\[PartialD]$$, $$0$$]1\)\), \ $$1000000000000000000$$])\)\)\)/(1 + \!$$\*SubscriptBox[\(\[PartialD]$$, $$0$$]$$\*SubscriptBox[\(\[PartialD]$$, $$0$$]1\)\)/1000000000000000000) (sorry I don't know how to show the formated output...) Jun 4, 2019 at 13:06
• I guess I see your point. I always used to do pot[y] /. y -> 0 // N , with a number as output. But the output I get from NDSolve does not change whether I use one or the other definition for the derivative, i.e., I still get the same constant plots... Jun 4, 2019 at 14:07