Multiple solutions using the shooting method

Question

$$\begin{cases}&-z^{\prime\prime}(t)=\lambda(1+(N-2)t)^{\frac1{2-N}(2(N-1)+\alpha)}f(z(t)),\quad t\in(0,+\infty)\\&z(0)=z^\prime(+\infty)=0\end{cases}$$

I'm trying to solve the above differential equation. I'm using the shooting method, but for each initial guess a different solution appears. Which solution is correct? Is it possible to monitor the convergence rate and the iteration steps of the shooting method? Any help or hint is welcome!

Here is the code:

NN = 3.;
p = 3.;
xf = 10000.;

sols = 
  Map[{z -> 
      NDSolveValue[{-z''[
           t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^
           p, z[0] == 0, z'[xf] == 0}, z, {t, 0, xf}, 
       Method -> {"Shooting", 
         "StartingInitialConditions" -> {z[0] == 0, z'[0] == #}}]} &,
   {0, 5, 10, 20, 40, 60, 140}];

Plot[Evaluate[z[t] /. sols], {t, 0, 100}, PlotRange -> All]

Answer 1

To assess the correctness of solutions, let's look at how well they meet the boundary conditions you impose.

derivativeAtZero = {0, 5, 10, 20, 40, 60, 140};

With[{NN = 3., p = 3., xf = 10000.}, 
  sols = 
      NDSolveValue[
        {-z''[t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^p,
         z[0] == 0, z'[xf] == 0},
        z, {t, 0, xf},
        Method ->
          {"Shooting",
           "StartingInitialConditions" -> {z[0] == 0, z'[0] == #}}] & 
      /@ derivativeAtZero];

computedDerivative = #'[0] & /@ sols;

TableForm[
  Transpose[
    {derivativeAtZero, computedDerivative, derivativeAtZero - computedDerivative}], 
  TableHeadings -> {None, {"Imposed", "Solved", "Difference"}}]

The 1st solution, which is the trivial solution z[t] = 0, looks good. For the rest, it's up to you decide if the interpolation function returned by NDSolveValue approximates the imposed boundary condition well-enough to satisfy you. The 2nd, 4th and 5th solutions don't look too bad, so let's plot them.

Plot[Evaluate[#[t] & /@ sols[[{2, 4, 5}]]], {t, 0, 20}, 
  PlotRange -> All, PlotLegends -> {2, 4, 5}]

On general grounds, only the plot corresponding to the 2nd boundary condition looks well behaved, but I can't be sure. I don't know the system being modeled, so I don't really know what kind of transient behavior can be expected.

Answer 2

There is no reason to suppose that the nonlinear ODE presented in the question has a unique solution. Indeed, three were identified here. Several more can be found with

NN = 3; p = 3; xf = 100000; 
sols = Map[{z -> 
    NDSolveValue[{-z''[t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^p, 
    z[0] == 0, z'[xf] == 0}, z, {t, 0, xf}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {z[0] == 0, z'[0] == #}}]} &,
   Range[0, 200, 5]];

and then identify the unique solutions.

DeleteDuplicates[sols, Abs[(z'[0] /. #1) - (z'[0] /. #2)] < .01 &];
Select[%, Evaluate[z'[0] /. #] > 0 &];
Sort[%, Evaluate[z'[0] /. #1] < Evaluate[z'[0] /. #2] &];
z'[0] /. %

(* {5.2103, 19.1366, 41.7079, 72.9223, 112.779, 161.279, 218.422, 284.207} *)

LogLinearPlot[Evaluate[z[t] /. %%], {t, 0.01, 1000}, PlotRange -> All]

Solutions have progressively more oscillations before becoming constant. In effect, we have a set of eigenfunctions.

It may be of interest to look at initial-value solutions for a Range of z'[0] near one of these eigenfunctions, for instance,

solt = Map[{z -> 
    NDSolveValue[{-z''[t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^p, 
    z[0] == 0, z'[0] == #}, z, {t, 0, 100000}]}&, Range[41.7075, 41.7085, .0001]];
LogLinearPlot[Evaluate[z[t] /. solt], {t, 0.01, xf}, PlotRange -> All]

demonstrating the sensitivity of the asymptotic solution to z'[0].