It is very easy to show (and everyone knows) that Nonlinear Schrodinger in the simplest case has Tanh kink solution:

f[x_] := Tanh[x/Sqrt[2]];
u''[x] + u[x] - u[x]^3 /. u -> f // FullSimplify

Out[] = 0

That is your equation at asymptotic 'x -> + Infinity' because

Limit[u'[x]/x - u[x]/x^2,x->0] == 0

And it is clear that solution around zero can be linearized. Now that's being said Micahel Trott had offered a proper way of finding this solution that you can see below in his reply to tis comment.

Yes, that sounds good. But what if you integtate to, say, x = 10? At least with my system the solution "explodes".

And it shows a peculiar behaviour around 1.11. Just have a try with 1.105 and 1.108 and so on....

Luca, Nice, you avoid the stiffness issue by using a shooting method! How did you find the initial value of 1.11?

Vitaly, Ok here it is:

In[201]:= ode1[u_] := D[u, {x, 2}] + D[u, x]/x - u/x^2 + u - u^3;

In[202]:= pdc[ode1[u[x]]]

Out[202]//TraditionalForm= -(u(x)/x^2)+\[PartialD]^2u(x)/\[PartialD]x^2-u(x)^3+u(x)+\[PartialD]u(x)/\[PartialD]x/x

transform equation

In[203]:= teq1 = Expand[ode1[u[x/(1 + x)]]] /. x/(1 + x) -> t /. x -> t/(1 - t)

Out[203]= u[t] - ((1 - t)^2 u[t])/t^2 - u[t]^3 + (
 2 t Derivative[1][u][t])/((1 - t) (1 + t/(1 - t))^3) - (
 3 Derivative[1][u][t])/(1 + t/(1 - t))^2 + ((1 - t) Derivative[1][u][t])/(
 t (1 + t/(1 - t))) + (
 t^2 (u^\[Prime]\[Prime])[t])/((1 - t)^2 (1 + t/(1 - t))^4) - (
 2 t (u^\[Prime]\[Prime])[t])/((1 - t) (1 + t/(1 - t))^3) + (
  u^\[Prime]\[Prime])[t]/(1 + t/(1 - t))^2

In[204]:= teq2 = Simplify[%] t^2

Out[204]= (-1 + 2 t) u[t] - 
 t^2 u[t]^3 + (-1 + 
    t)^3 t ((-1 + 2 t) Derivative[1][u][t] + (-1 + t) t (u^\[Prime]\[Prime])[
      t])

In[205]:= Factor[teq2]

Out[205]= -u[t] + 2 t u[t] - t^2 u[t]^3 + t Derivative[1][u][t] - 
 5 t^2 Derivative[1][u][t] + 9 t^3 Derivative[1][u][t] - 
 7 t^4 Derivative[1][u][t] + 2 t^5 Derivative[1][u][t] + 
 t^2 (u^\[Prime]\[Prime])[t] - 4 t^3 (u^\[Prime]\[Prime])[t] + 
 6 t^4 (u^\[Prime]\[Prime])[t] - 4 t^5 (u^\[Prime]\[Prime])[t] + 
 t^6 (u^\[Prime]\[Prime])[t]

In[206]:= eps = .1;

In[207]:= NDSolve[{teq2 == 0, u[0] == 0, u[1.] == 1.}, u[t], {t, eps, 1.}]

During evaluation of In[207]:= Power::infy: Infinite expression 1/0.^2 encountered. >>

During evaluation of In[207]:= Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

During evaluation of In[207]:= Power::infy: Infinite expression 1/0.^2 encountered. >>

During evaluation of In[207]:= Power::infy: Infinite expression 1/0.^2 encountered. >>

During evaluation of In[207]:= General::stop: Further output of Power::infy will be suppressed during this calculation. >>

During evaluation of In[207]:= Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

During evaluation of In[207]:= Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

During evaluation of In[207]:= General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

During evaluation of In[207]:= NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

Out[207]= NDSolve[{(-1 + 2 t) u[t] - 
    t^2 u[t]^3 + (-1 + 
       t)^3 t ((-1 + 2 t) Derivative[1][u][t] + (-1 + t) t (
         u^\[Prime]\[Prime])[t]) == 0, u[0] == 0, u[1.] == 1.}, 
 u[t], {t, 0.1, 1.}]

This blows up because of division by zero at the low end. I can fix that by using the low end asymptotic term as t->0 t u'==u:

In[217]:= NDSolve[{teq2 == 0, t u'[eps] == u[eps], u[1.] == 1.}, 
 u[t], {t, eps, 1.}]

During evaluation of In[217]:= NDSolve::ndsz: At t == 0.9999112864801112`, step size is effectively zero; singularity or stiff system suspected. >>

Out[217]= NDSolve[{(-1 + 2 t) u[t] - 
    t^2 u[t]^3 + (-1 + 
       t)^3 t ((-1 + 2 t) Derivative[1][u][t] + (-1 + t) t (
         u^\[Prime]\[Prime])[t]) == 0, 
  t Derivative[1][u][0.1] == u[0.1], u[1.] == 1.}, u[t], {t, 0.1, 1.}]

Now it blows up near t=1 because the differential equation turns into an algebraic equation (u=u^3) via an apparent stiffness at 1 (coefficients of the derivatives start to vanish)

Attachments:

Actually, Kay, I like your transform much better. Could you post your code or attach a notebook or the transformed ODE and NDSolve errors? We could try to figure out a custom method in NDSolve that would handle stiffness.

Vitaly, Your transform is similar to the one I tried t=x/(1+x) yours: t=Sqrt[x]/(1+Sqrt[x]) . Any reason Sqrt[x}? In any case, NDSolve was unable to solve the transformed equation. I think the reason is that the equation becomes very stiff near the boundaries. Although I'm sure there is a way, I didn't really spend a lot of time on it.
Mattia mentioned that the equation was solved before, maybe it is worth looking at some literature to see how the problem has been handled in the past.

Ah, ok! When deriving the new domain, I had a typo and used $t = e^{-x}$ instead of $x = e^{-t}$...

I attacked this problem by making the substitution $x = e^{-t}$. When subbing for the independent variable, we need to make careful use of the chain rule to express $u(x), u'(x), u''(x)$ in terms of $u(t), u'(t), u''(t)$.

ode = u''[x]+(u'[x]/x)-(u[x]/(x^2))+u[x]-u[x]^3 == 0;
sub = Exp[-t];

newode = ode /. {
    u[x] -> u[t],
    u'[x] -> u'[t]/D[sub, t],
    u''[x] -> (u''[t] - D[sub, {t, 2}] u'[t]/D[sub, t])/D[sub, t]^2,
    x -> sub
}

u[t] - E^(2 t) u[t] - u[t]^3 - E^(2 t) u'[t] + E^(2 t) (u'[t] + u''[t]) == 0

sol = NDSolveValue[{newode, u[1] == 0, u[0] == 1}, u[t], t] /. t -> Exp[-x];
sol /. x -> 0

0.

sol /. x -> 10^15

1.

Plot[Evaluate[sol], {x, 0, 10}, PlotRange -> All]

Here's my issue though... If I plug this solution into the original ODE, I don't get zero, so I think I did something wrong but I can't figure out where. In fact I am having the same issue with Hans Dolhaine's solution he posted above (his is sol2 in the plot below).

Maybe someone here can figure out where I went wrong? I feel it's a simple oversight.

Plot[Evaluate[{
   D[sol, {x, 2}] + (D[sol, x]/x) - (sol/x^2) + sol - sol^3,
   D[sol2, {x, 2}] + (D[sol2, x]/x) - (sol2/x^2) + sol2 - sol2^3
 }], {x, 4, 20}, PlotRange -> All]

You do not have to repeat the problem to me, I've read it. The idea that I am suggesting is that the statement

"I don't think that domain has necessary to be infinitive but big enough"

leads to an approximation of BC at Infinity - nothing wrong with that. I am just saying potentially there is a mapping of domain that would allow numerical routine to run on the full interval.

NDSolve does not take Infinity as a spec for BC. So one has to trick out, use asymptotic, approximations. Nothing wrong with that - I am just suggesting a potential way to avoid that and rely completely on numerical routine via easy BC obtained from domain mapping.

This mapping could be wrong too because of potential issues that could appear during solution. Maybe other mapping can be tried then. Just an idea.

Folks I have not check this to the end but I have a question to all of you. One of the tight spots here is the BC at Infinity. I've seen someone mention mapping (0, Infinity) -> (0, 1). Have you looked into it? There is a beautiful algebraic mapping for the positive half-plane to (0,1):

t[x_] := 1 /(1 + x^(-1/2));
Plot[t[x], {x, 0, 10}, PlotRange -> All]

Series[t[x], {x, 0, 1}]
Series[t[x], {x, Infinity, 1}]

Has anyone tried this change of variable to be at ease using BC now at (0,1) using NDSolve ?

The new equation is:

ode1[u_] := 
  D[u, {x, 2}] + D[u, x]/x - u/x^2 + u - u^3 - 
    D[x*D[u^2, x], x]*u/x  == 0;

After the transformation and considering only dominating terms, I tried to use:

DSolve[u1''[t] - u1[t] - u1[t]*(u1'[t])^2 - 2*u1''[t]*(u1[t])^2 == 0, 
 u1[t], t]

but I had no output...

keep in mind, my approximation is very rough and not accurate. I'm sure there is a way to get a better solution

Dear Hans, it's a very peculiar situation. But it is possible to show that for a large system with a fix number of particles, the state of lowest energy is a uniform distribution. In this calculation, the effect of the boundaries may be neglected. You can consider the system as "unbounded". This is a n-body wave function, for n identical particles (an imperfect Bose gas). Ginzburg and Pitaevskii solved numerically that equation in 1958, I don't know how, and showed that the solution is a monotone curve...

The other equation is the same with another term and a parameter: $$u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3 - b*(1/x)*D[x*D[u[x]^2]*u[x]=0$$. I hope to be able to get a numerical solution varying b using your idea.

Thanks a lot Kay!

Dear Hans, unfortunately the behavior of the first derivative is unknown. This is a typical Schroedinger problem, in which only BC's are known. Your solution seems to be correct near 0, but the function is monotone and most reach asymptotically 1 at +inf.

Hello Mattia,

as it is a diff.eq. of second order we need to have information about the first derivative as well. And the division by x leads to problems at x == 0. But we could try to start in the vicinity of 0:

eq = u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3 == 0;
eps = 10^-4;
sol = u /. 
  NDSolve[{eq, u[eps] == 0, u'[eps] == .2}, u, {x, eps, 5}][[1]]
Plot[sol[x], {x, eps, 3}]

If you assume that the solution is u = x^(Sqrt[5]/2 - 1/2)/(B + x^(Sqrt[5]/2 - 1/2))

and the equation is exact at x=1 then B -> (-3 + Sqrt[5] - Sqrt[38 - 14 Sqrt[5]])/(2 (-3 + Sqrt[5])) or about 2.19353 or approximation u=x^0.618/(2.1935+x^0.618)

which satisfies your boundary conditions