Imposing decaying boundary conditions on a non linear ode

I am trying to solve $$a^{2}y''=y+y^3$$ numerically. This equation models a potential and goes to $$\infty$$ for $$x\to0$$ hence I get the singularity to be of order $$\frac{1}{x}$$ by keeping only the y^{3} term similarly decays exponentially which can be seen by keeping only the y term. I thought that it would be necessary to remove the singularity for any algorithm to be stable. So I converted the Ode using the substitution u=x^3 y which I used to make sure it is a also differentiable at zero. What I can't figure out is how to actually implement the boundary conditions numerically which are that it vanish at $$\pm \infty$$ and the singularity is fixed to be at the origin or a zero in the case of the transformed ode. What boundary value should I actually use Numerically? If I fix both boundary points to have zero value then it gives the trivial solution zero? Also it would be great if someone could tell me whether there are numerical algorithms for dealing with odes with singularities?

'

• From what you described it is impossible to say why exactly you don't see a mesh convergence by decreasing step size... Please put some code as well... – Alone Programmer Oct 31 '19 at 14:44
• @Aloneprogrammer 8 have added the code it's in python. I think the problem here is that I am giving the wrong initial conditions I use close to zero but I am not sure what that means 0. 1 something smaller. The main problem is that it is technically a boundary value problem but x=0 singularity interferes with those algorithms – ben tenyson Nov 2 '19 at 5:32
• 1. How can you solve an ODE without specifying initial or boundary conditions? 2. What goes to infinity as $x\to0$? The ODE as written is autonomous so there is nothing special about $x=0$. – Rahul Nov 2 '19 at 11:35
• @Rahul I understand that the Ode is autonomous which implies that I can shift the singularity to any point.The exact solution for the case of x=0 pole and a=1 is $\sqrt{2}\frac{1}{\sinh(x)}$. I want to numerically obtain this result. I transformed the Ode using substitution described above to remove the singularity from the solution. After doing, I do get the value of the y and y' at x=0 which is 0 because of the way I transformed it but the new ode is not continuous at x=0 so that doesn't work then I have started to trying different initial value away from x=0 and I get the right shape – ben tenyson Nov 2 '19 at 12:04
• @Rahul but there is a consistent difference between the exact and numerical solution – ben tenyson Nov 2 '19 at 12:05