I am trying to solve the non linear ODE
$$\displaystyle \frac{d^2y} {dx^2}=a^2(y+y^3)$$
With the boundary conditions that it vanish at $\pm\infty$. I am thinking that it might be better to deal with it an integral equation
$$\displaystyle \left(\frac{d^2} {dx^2}-a^2 \right) y =a^2y^3$$
The Green's function will be for that of the modified Hemholtz equation in 1D and the after stacking the solution to the homogenous equation(following approach of Lippmann Schwinger equation) I can write it in integral form as
$$\displaystyle y(x)=e^{-a|x|}-\frac{1}{2a}\int_{-\infty}^{\infty}e^{-a|x-x_{1}|}y^{3}(x_{1})dx_{1}$$
Now I have no idea how to solve this integral equation. For the Lippmann Schwinger equation Neumann series can be used should I do the same with this one or is there some better method for these kind of integral equations?
You can integrate once, after multiplying with $2y'$ $$ y'^2=a^2y^2+\frac{a^2}2y^4+c $$ Vanishing at infinity implies coming to rest there, so that the solution that you seek has $c=0$. Then $$ y'=\pm ay\sqrt{1+\frac12y^2} $$ This has the trivial solution $y=0$. No other solution in both sign variants can change its sign. All further considerations assume that the solution is either entirely positive or entirely negative.
Substitute $y=\sqrt2\sinh(u)$ and assume $y\ne 0$, $u\ne 0$ $$ \sqrt2\cosh(u)u'=\pm a\sqrt2\sinh(u)\cosh(u)\implies \frac{2e^uu'}{e^{2u}-1}=\pm a \\~\\ \frac{e^u-1}{e^u+1}=Ce^{\pm ax}\implies e^u=\frac{1+Ce^{\pm ax}}{1-Ce^{\pm ax}} $$ so that finally $$ y=\frac{e^u-e^{-u}}{\sqrt2}=\frac{2\sqrt2Ce^{\pm ax}}{1-C^2e^{\pm 2ax}}=\pm\frac{\sqrt2}{\sinh(a(x-x_0))} $$
One could restart now with $y(x)=\frac{\sqrt{2}}{\sinh(v(x))}$ to get $$ \frac{2\cosh^2(v)v'^2}{\sinh^4(v)}=\frac{2a^2}{\sinh^2(v)}\left(1+\frac1{\sinh^2(v)}\right)\implies v'^2=a^2, $$ which indeed directly gives $v(x)=\pm a(x-x_0)$.
As any of these solutions has a pole at $x=x_0$, they do not solve the given problem, the stationary solution $y=0$ is the only one satisfying all conditions.