I want to solve a differential equation numerically and determine one of the parameters in the D.E. The goal is to determine $M$ such that the function goes asymptotically to zero or to put it in a better way close to zero, without ever getting negative values.
What I want to solve is the following
$$f''(r) + \frac{6}{r} f'(r) - \frac{1}{r^2} \Big(\frac{1}{2}\Big)^2 f(r) + \frac{6}{r^2}f(r) + \frac{1}{r^2} \frac{1}{2} f(r) + \frac{M^2}{(r^2 + lf(r)^2)^2}f(r) = 0 $$
In order to calculate the $lf(r)$ in the above, some calculations are needed. Below, I give the code and I am breaking it in two pieces. The first is what is needed to calculate $lf(r)$ and the second is my attempt to solve the D.E.
The first piece of the code. I don't have any questions here and this is just to evaluate $lf(r)$ in the above D.E.
rmax := 10^4;
nc := 2;
nf1 := 2;
nf2 := 4;
\[Alpha]1 := 0.7;
b0 := 1/(6 \[Pi]) (11 (nc + 1) - 2 nf1 - 2 nf2 (nc - 1));
b1 := 1/(24 \[Pi]^2) (34 (nc + 1)^2 - (10 (nc + 1) +
6*(2 nc + 1)/4) 1/2 2 nf1 - (10 (nc + 1) + 6 nc) (nc - 1) nf2);
b0;
b1;
eqalp1 = alp1'[mu1] + (b0 alp1[mu1]^2 + b1 alp1[mu1]^3) == 0;
eqalp1;
sol1 = NDSolve [{ eqalp1, alp1[0] == \[Alpha]1},
alp1, {mu1, 0, 10^24}];
alph1[mu1_] := First[Evaluate[alp1[mu1] /. sol1]]
alph1[mu1];
\[Alpha][mu_] := alph1[mu]
\[Gamma][mu1_] := 3/(2 \[Pi]) \[Alpha][mu1] ((2 nc + 1)/4 + nc)
\[Delta][mu1_] := -2 \[Gamma][mu1]
muv = 0;
ode := D[r^3 L'[r], r] - \[Delta][Log[Sqrt[r^2 + L[r]^2] ]] r L[r] == 0
solm[v_] := NDSolve[{ode, L[v ] == v,
L'[v] == 0}, L, {r, v, rmax}, Method -> {"StiffnessSwitching"},
MaxSteps -> Infinity, PrecisionGoal -> 13, AccuracyGoal -> 13] //
First
LinfBm[ v_?NumericQ ] := (L[rmax] /. solm[v]) - muv
mir = m /. FindRoot[LinfBm[m], {m, 1}]
lf[\[Rho]_] := Evaluate[L[\[Rho]] /. solm[ mir]]
And the second part, the one that performs the NDSolve.
diracpde =
D[f[r], {r, 2}] + 6/r D[f[r], r] - 1/r^2 (1/2)^2 f[r] +
6 /r^2 f[r] + 1/r^2 (1/2) f[r] + M^2/(r^2 + lf[r]^2)^2 f[r];
diracforndsolve = {diracpde == 0, f[mir] == 1, f'[mir] == 0};
sol1 = NDSolve[diracforndsolve /. M :> 2, f[r], {r, mir, rmax}];
Plot[f[r] /. sol1, {r, mir, rmax},
PlotRange -> {{mir, rmax}, {-1, 1}},
BaseStyle -> {18, FontFamily -> "Times New Roman"},
AxesLabel -> {"\[Rho]", "top(\[Rho])"}, PlotStyle -> {Thick, Red}]
As mentioned above, the goal is to change the value of $M$ and find the solution that goes to zero, but never crosses to negative values.
One way of doing so, but I don't really like it is to plot from zero to one instead of -1 to 1.
My question is, if there is a better way to test it.
I have tried the following
h[r_] := f[r] /. sol1
FullSimplify[ForAll[r, mir <= r <= rmax, h[r] > 0]]
h[r_] := f[r] /. sol1
Reduce[ForAll[r, mir <= r <= rmax, h[r] > 0]]
FullSimplify[ForAll[r, mir <= r <= rmax \[Implies] h[r] > 0]]
motivated by some answers I found in other threads, but they don't do the job.
Another question I have been having is the following. If I set a lower $rmax$, i.e if I set it to $10$, the D.E has a solution with the feature tht I want and can be easily found, but for the case above $10^4$ it seems that it does not, or that I have done something wrong.
Any help/comments are welcome.