利用ParametricNDSolveValue求解非线性微分方程组

我在利用ParametricNDSolveValue求解非线性微分方程组时得不到我想要解，“h''[r] + (3 \[Rho]'[r]/\[Rho][r] h'[r] - dV[h[r]]) == 0, \[Rho]''[r] + k1*\[Rho][r]/3 (h'[r]^2 + V[h[r]]) == 0”这是我的方程组，其中“k1 = 1; V[h_] := 235.17737975645832 h^2 - 0.028362365242459223 h^4 + h^6/2976800; dV[h_] = D[V[h], h];”然后它的边界条件是“h'[10^-2] == 10^-2, \[Rho][10^-2] == 10^-2, \[Rho]'[10^-2] == 1, 0 <h[rmax] < 10^-3”这里的rmax是一个有限的距离，在rmax处h[r]无限趋于0。我以下的程序采用的是打靶法来解决这个问题。但他并不能实现h[r]在rmax处无限趋于0。但当我取”V[h_] := 0.1 + 0.175 h^2 - 0.45 h^3 + h^4/4; h0Start = 0.80621; h0Step = 0.00001; maxAttempts = 10;“时，他画出来的图可以得到满足我的需要的解。
程序如下：
$Assumptions = h > -1;
V[h_] :=
235.17737975645832` h^2 - 0.028362365242459223` h^4 + h^6/2976800;
dV[h_] = D[V[h], h]; rmin = 10^-3; k1 = 1;
sys = {h''[r] + (3 \[Rho]'[r]/\[Rho][r] h'[r] - dV[h[r]]) ==
0, \[Rho]''[r] + (k1/3) \[Rho][r] (h'[r]^2 + V[h[r]]) == 0,
h[rmin] == h0,
h'[rmin] == 10^-2, \[Rho][rmin] == 10^-2, \[Rho]'[rmin] == 1,
WhenEvent[h[r] < 10^-30, "StopIntegration"]};
pf = ParametricNDSolveValue[sys, {h, \[Rho]}, {r, rmin, 500}, {h0},
Method -> {"StiffnessSwitching"}, MaxSteps -> 10^6,
PrecisionGoal -> 10, AccuracyGoal -> 10, "MaxStepSize" -> 0.001]
CheckSolutionQuality[h0_?NumericQ] :=
Module[{sol, hsol, rmaxActual, testPoints, derivatives},
sol = Quiet@pf[h0];
hsol = sol[[1]];
rmaxActual = hsol["Domain"][[1, 2]];
If[rmaxActual >= 500 - 1,
Return[{"IntegrationFailed", h0, hsol, sol[[2]]}]];
If[hsol[rmaxActual] > 10^-3,
Return[{"FailAtRmax", h0, hsol, sol[[2]]}]];
testPoints = Subdivide[RH + 0.01, rmaxActual, 200];
derivatives = hsol' /@ testPoints;
Which[AllTrue[derivatives, # < -10^-20 &],
Return[{"StrictlyDecreasing", h0, hsol, sol[[2]]}],
AllTrue[derivatives, # <= 0 &],
Return[{"NonStrictDecreasing", h0, hsol, sol[[2]]}], True,
Return[{"NotMonotonic", h0, hsol, sol[[2]]}]]]
h0Start = 80;(*找适合边界条件h[rmax]趋于0的初始值，因为按物理情景来说他应该比较靠近h0=227这个方向*)
h0Step = 0.00001;
maxAttempts = 10;
strictMode = False;
h0failed = h0Start;
result = Catch[Do[currentH0 = h0Start + k*h0Step;
status = CheckSolutionQuality[currentH0];
(*记录最后一次失败的 h0*)
If[status[[1]] =!= "StrictlyDecreasing" &&
status[[1]] =!= "NonStrictDecreasing", h0failed = currentH0;];
Which[
status[[1]] ==
"StrictlyDecreasing" || (status[[1]] ==
"NonStrictDecreasing" && ! strictMode),
Throw[{"Success", currentH0, status[[1]]}],
status[[1]] == "IntegrationFailed",
Print["警告: h0=", currentH0, " 积分失败，尝试调整步长"];, True,
Print["尝试 h0=", currentH0, " 失败: ", status[[1]]];
solFailed = Quiet@pf[currentH0];
hSolFailed = solFailed[[1]];
rmaxFailed = hSolFailed["Domain"][[1, 2]];
If[rmaxFailed + 0.1 <= 100,
Print[" rmax 的值: ", rmaxFailed, " 解在 rmax 处的值: ",
hSolFailed[rmaxFailed], "解在 rmax+0.1 处的值: ",
hSolFailed[rmaxFailed + 0.1]],
Print["警告: rmax+0.1 超出范围，无法计算"]];], {k, 0, maxAttempts}];
{"Failed", "超过最大尝试次数"}];
(*结果处理与可视化*)
If[result[[1]] === "Final Success" || result[[1]] === "Success",
h0opt = result[[2]];
{hSolution, \[Rho]Solution} = pf[h0opt];
rmax = hSolution["Domain"][[1, 2]];
Print["========================================"];
Print["成功找到解！"];
Print["- 最佳 h0 = ", h0opt];
Print["- rmax = ", rmax];
Print["- 最终 h(rmax) = ", hSolution[rmax]];
Print["- 单调性模式: ", result[[3]]];
Print["- 最小h值: ", NMinimize[{hSolution[r], 10^-3 <= r <= rmax}, r]];
(*绘制成功的解*)
Grid[{{Plot[hSolution[r], {r, RH + 10^-6, rmax},
PlotLabel -> "场分布 h(r)", ImageSize -> 300, PlotRange -> All],
Plot[hSolution'[r], {r, RH + 10^-6, rmax},
PlotLabel -> "导数 h'(r)", ImageSize -> 300,
PlotRange -> All]}, {Plot[\[Rho]Solution[r], {r, RH + 10^-6,
rmax}, PlotLabel -> "质量分布 \[Rho](r)", ImageSize -> 300,
PlotRange -> All], SpanFromLeft}}, Frame -> All],
Print["========================================"];
Print["未能在 ", maxAttempts, " 次尝试内找到有效解"];
Print["最后尝试的 h0 = ", h0failed];
Print["失败原因: ", result[[2]]];
Print["解在rmax+1处的值: ", hSolFailed[rmaxFailed + 0.1]];
solFailed = Quiet@pf[h0failed];
hSolFailed = Quiet@pf[h0failed][[1]];
\[Rho]SolFailed = solFailed[[2]];
rmaxFailed = hSolFailed["Domain"][[1, 2]];
(*绘制失败的解*)
Grid[{{Plot[hSolFailed[r], {r, RH + 10^-6, rmaxFailed},
PlotLabel -> "失败解: 场分布 h(r)", ImageSize -> 300,
PlotRange -> All],
Plot[hSolFailed'[r], {r, RH + 10^-6, rmaxFailed},
PlotLabel -> "失败解: 导数 h'(r)", ImageSize -> 300,
PlotRange -> All]},
{Plot[\[Rho]SolFailed[r], {r, RH + 10^-6, rmaxFailed},
PlotLabel -> "失败解: \[Rho](r)", ImageSize -> 300,
PlotRange -> All], SpanFromLeft}}, Frame -> All]]