Non-unique solutions for ODE with boundary conditions at infinity Ello,
I am looking for solutions to equations such as
$$ U ^ \prime (y) - y^{\prime \prime} + \beta y^{\textit{IV}} = 0$$
where $\beta >0$ is a constant and $U(x)$ is a function whose Taylor expansion at the origin equals $ \alpha x^2 + \dots$ for $\alpha >0$ and such that $\lim_{x \to \infty} = -\infty$.
A prototypical example would be the quartic double-well, $$ U(x) = x^2 -x^4$$.
There is a trivial solution $y(x) = 0$. On intuitive grounds (please see physical background below), as well as noting that the origin is a hyperbolic point for a dynamical system governed by the related system of equations, another solution is expected as well, with boundary conditions $$ \begin{cases} \lim_{x \pm \infty} y(x) = 0  \\ 
\lim_{x \pm \infty} y^ \prime (x) = 0 \end{cases}$$
This fact does not contradict the unicity theorem, as far as I understand, as the latter deals with boundary conditions on the real line, I believe.
Now to my problem, what numerical scheme could be used to get the non-trivial solution? Anything I could think of will pick the trivial solution branch only.
Further, I would be grateful for theory references on the possibility of multiple solutions for boundary conditions at infinity.
Things would be different if $\beta$ were equal to $0$.
In that case it is easy to determine two boundary conditions (also clarified in the physical background section below) for $x = 0$, and any numerical method will handle it easily.
Physical Background
The ODE is the Euler-Lagrange equation for the functional
$$ \int _{- \infty} ^{\infty } \frac{1}{2} y^{\prime 2} + U(y) + \frac{\beta}{2} y^{\prime \prime 2} \mathrm{d}x $$
representing the energy of a stretchable beam in the potential $U$. The first term represents the stretching energy, the second the potential energy, the third the bending energy.
Given this physical picture, the trivial solution represents a beam lying in the potential energy "valley" at the origin, un-stretched and straight. The non-trivial solution represents a beam lying on the valley as before, then deforming over the potential energy "bump", kept in the equilibrium by the pull from the descending potential energy on the other side, and going back to the valley at the origin.
If $\beta = 0$ the non-trivial solution is amenable to a simple solution.
The energy functional reduces in such case to
$$ \int _{- \infty} ^{\infty } \frac{1}{2} y^{\prime 2} + U(y)  \mathrm{d}x $$ which could be interpreted as the action of a particle rolling in a potential energy $-U(x)$ (with $x$ now representing time). The turning point is then easily found by solving the equation $ U(x_1) = 0$: that is, the particle starts from $x = 0$, rolls down and up until $x_1$, and rolls back in a homoclinic orbit.
This gives two boundary conditions
$$ \begin{cases}  y(0) = x_1  \\ 
 y^ \prime (0) = 0 \end{cases}$$
and these allow any numerical scheme to pick the non-trivial branch.
If $\beta \neq 0$, this simple reasoning fails for the full fourth-order ODE (I also thought about converting the fourth-order ODE to a system of equations, to attempt a similar interpretation for a system of four particles, but did not get anywhere yet, apart confirming that a second non-trivial solution exists). I am hence left with the boundary conditions a infinity, and cannot get any solution other than the trivial one.
To sum up, I would like to be pointed towards references for the non-unicity of solutions when boundary conditions at infinity are applied, and to understand how to numerically solve similar problems.
The ODE has a physical interpretation and I was unsure whether the question were more apt for Physics StackExchange, or even SciComp StackExchange. I would certainly re-route if advised, I thought to start here as ultimately the question is primarily on the maths, thanks a lot.
 A: In the far field at $y=+\infty$ you have that everything is zero, so on the approach you get small values. This allows to ignore the higher-degree terms in $U(y)$ resp. $U'(y)$ and reduce the equation to an approximate linear DE
$$
βy^{(4)}-y''+2αy=0.
$$
This has a characteristic polynomial
$$
0=βλ^4-λ^2+2α\implies (βλ^2+\sqrt{2αβ})^2=β(1+2\sqrt{2αβ})λ^2
$$
etc. Better find the roots numerically.
Looking closely, one can read off that the solutions come in pairs with positive and negative real part. For the part to positive infinity the coefficients to the roots with positive real part have to be zero. This can be translated into a boundary condition at some large finite point.
The same can be done for the far field in direction of negative infinity.
So taking the 4 roots $\lambda_{\pm k}$, $k=1,2$, with $Re(λ_{+k})>0$ etc. then the far field pieces are for $x\ge T$
$$
y_+(x)=A_+e^{λ_{-1}x}+B_+e^{λ_{-2}x}
$$
and for $x<-T$
$$
y_-(x)=A_-e^{λ_{+1}x}+B_-e^{λ_{+2}x}
$$
Leaving the coefficients as free parameters, or as constant components in the ODE system, gives a total of 8 components and thus 8 slots for boundary conditions at $\pm T$.
With these one can establish the general boundary value problem with the first 6 boundary conditions being the trivial
\begin{align}
y(T)&=y_+(T),& y(-T)&=y_-(-T) \\
y'(T)&=y_+'(T),& y'(-T)&=y_-'(-T) \\
y''(T)&=y_+''(T),& y''(-T)&=y_-''(-T) \\
\end{align}
Using the remaining 2 slots to ensure 3rd order derivative continuity would force the solution to be the zero function. Instead one could use one slot to fix the integral over $y^2$ or similar and use the remaining degree of freedom to minimize the jumps in the 3rd derivative. This should give convincing solution approximations.

(added 2022/4/2) Continuing from the partial factorization, completing it gives factors
$$
βλ^2\pm \sqrt{β(1+2\sqrt{2αβ})}λ+\sqrt{2αβ}=0.
$$
Applying Viete rules tells that the version with the positive sign combines the roots with negative real part and the version with negative sign the roots with positive real part. The corresponding second order DE are those that the solution approximately follows going to $\pm\infty$, respectively, if it stays bounded there. Thus one can get sensible boundary conditions without introducing additional parameters as
$$
βy''(T)+ \sqrt{β(1+2\sqrt{2αβ})}y'(T)+\sqrt{2αβ}y(T)=0
\\
βy''(-T)- \sqrt{β(1+2\sqrt{2αβ})}y'(-T)+\sqrt{2αβ}y(-T)=0
$$
and their versions in one derivative higher.
To force a non-zero solution add the energy integral as
$$
I'(x)=Ky(x)^2,~~ I(-T)=-1,~~I(T)=1.
$$
The parameter $K$ is also a free parameter or a constant state component, depending on what the solver permits to use.
A: You can use a linearization procedure to solve this nonlinear boundary value. For instance using a linear approximation scheme consisting in iterations over the linear ODE
$$
y''_{k+1}-V(y_k)+V'(y_k)(y_{k+1}-y_k)-\beta y''''_{k+1}=0
$$
given $y_0$ and the boundary conditions
$$
\cases{
y_{k+1}(-x_{max})=0\\
y'_{k+1}(x_{max})=0\\
y''_{k+1}(-x_{max})=0\\
y''_{k+1}(x_{max})=0
}
$$
Let us proceed with the calculations for $\beta = 0$ which is known to have as solution
$$
y=\sqrt{1-\tanh ^2\left(\sqrt{2} x\right)}
$$
We will approach this solution using negative $\beta$'s as follows in the attached MATHEMATICA script.
Clear[V, dV, y1, y0]
V[x_] := -2 x + 4 x^3;
dV[x_] := -2 + 12 x^2

y0 = Exp[-x];
beta = -0.00001;
xmax = 5;
nmax = 20;
bcs = {y1[-xmax] == 0, y1[xmax] == 0, y1'[-xmax] == 0, y1'[xmax] == 0};
SOLS = {Plot[y0, {x, -xmax, xmax}]};
thickness = Thin;
color = Blue;

For[k = 1, k <= nmax, k++,
   ode = y1''[x] + V[y0] + dV[y0] (y1[x] - y0) - beta y1''''[x];
   sol = NDSolve[Join[{ode == 0}, bcs], y1, {x, -xmax, xmax}][[1]];
   y0 = y1[x] /. sol;
   If[k == nmax, thickness = Thick; color = Black];
   AppendTo[SOLS, Plot[y0, {x, -xmax, xmax}, PlotStyle -> {thickness, color}]]
]

grfin = SOLS[[Length[SOLS]]];
Show[Plot[Sqrt[1 - Tanh[Sqrt[2] x]^2], {x, -xmax, xmax}, PlotStyle -> {Thick, Red}], grfin]
Show[SOLS, grfin, PlotRange -> All]

so with $\beta = -0.00001$ we obtain

