How to detect the positive solution to this logistic boundary value problem

Question

I am trying to find the positive solution to this nonlinear ODE. Theoretically, it has been proven that this problem has a unique positive solution for lam large. But, Mathematica can detect only the zero solution (zero is a solution for all lam's). Any techniques to detect the positive solution?

lam = 30;
sol = NDSolve[{-D[f[x], x, x] == lam*f[x]*(1 - f[x]), 
(D[f[x], x] /. x -> 1) == -Sqrt[lam]*0.1*f[1], (*first boundary condition*) 
(D[f[x], x] /. x -> 0) == Sqrt[lam]*0.1*f[0]}, (*second boundary condition*) 
f, {x, 0, 1}];

Plot[f[x] /. sol, {x, 0, 1}, AxesLabel -> {x, f}, LabelStyle -> Directive[Black, Bold, 12]]

Answer 1

Let us try a shooting method. The following code produces the solution for a given value of f[0]:

lam = 30;
sol[f0_?NumericQ]:=NDSolveValue[{-f''[x]==lam*f[x]*(1-f[x]),
                                 f[0]==f0,f'[0]==Sqrt[lam]*0.1*f0},f,{x,0,1}];

We want to understand for which f[0] the boundary condition at x==1 is satisfied. This condition is that the following must be zero:

cond1[f_]:=f'[1]-(-Sqrt[lam]*0.1*f[1])

Let us plot this:

Plot[cond1[sol[f0]],{f0,-0.53,0.91},AxesLabel->{"f[0]","cond1"}]

I picked the plot domain after some playing around. We see four values of f[0] where the condition at x==1 is also satisfied. I think the one OP asks for is the one on the very right, so let us see where it is:

f0x=f0/.FindRoot[cond1[sol[f0]],{f0,0.9}]
(* 0.905557 *)

Plot the corresponding solution:

With[{s=sol[f0x]}, Plot[s[x],{x,0,1},AxesLabel->{"x","f[x]"}]]

It is everywhere positive. It is also not difficult to guess that this solution is symmetric, $f(x) = f(1-x)$, and a numerical check confirms this. No error or warning message was produced, so I am assuming we can trust the numerical evaluation. But of course, one should check.

Note. The above is a manual implementation of the shooting method. Actually NDSolve has a shooting method built in. To see how this works, see this tutorial. Here:

lam=30;
sol2[f0guess_?NumericQ]:=NDSolveValue[{
                     -f''[x]==lam*f[x]*(1-f[x]),
                     f'[0]==Sqrt[lam]*0.1*f[0],
                     f'[1]==-Sqrt[lam]*0.1*f[1]},
                     f,x,Method->{"Shooting","StartingInitialConditions"->
                          {f[0]==f0guess,f'[0]==Sqrt[lam]*0.1*f0guess}}];

Let us use this to plot the four solutions mentioned above:

With[{sols=Map[sol2[#][x]&,{-0.5,0,0.85,0.9}]},
  Plot[Evaluate[sols],{x,0,1},PlotStyle->{Black,Blue,Green,Magenta},AxesLabel->{"x","f[x]"}]]

Note 2. The quantity 1/2*f'[x]^2+lam*(1/2*f[x]^2-1/3*f[x]^3) is independent of x for every solution.