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

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, (*first boundary condition*)
(D[f[x], x] /. x -> 0) == Sqrt[lam]*0.1*f}, (*second boundary condition*)
f, {x, 0, 1}];

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


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

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


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

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


Let us plot this:

Plot[cond1[sol[f0]],{f0,-0.53,0.91},AxesLabel->{"f","cond1"}] I picked the plot domain after some playing around. We see four values of f 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'==Sqrt[lam]*0.1*f,
f'==-Sqrt[lam]*0.1*f},
f,x,Method->{"Shooting","StartingInitialConditions"->
{f==f0guess,f'==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.