# FindMaximum not working with ParametricNDSolve

I have the following code (ignore the manipulate parameters, assume them constant in this context).

For background, what I'm trying to do is take the results of a 2D differential equation, find where it intersects y=0, and maximize the x where that happens for different parameter values.

If I replace FindMaximum with Plot (and change {R, 0.05} to {R, 0.005, 0.05}) I get a perfect plot of the x intercept as a function of R, with a visible maximum around R=0.02. However, the below code throws a lot of errors, including

FindRoot::nlnum: The function value {ParametricFunction[................][R]} is not a list of numbers with dimensions {1} at {t} = {40.}.
ReplaceAll::reps: {FindRoot[(y[t]/. sol$40848)[R],{t,40,0.05,50}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. FindMaximum::nrnum: The function value -<<1>> is not a real number at {R} = {0.05}.  In some iterations of this, when I was using NDSolve, I was actually getting in an infinite loop where the kernel would hang on spewing repetitive error messages (mostly about R not being replaced in NDSolve). What's going on? ρm = 1600 ;(*Density of rock*) ρa = 1.225;(*Density of air*) Cd = 0.5;(*Typical drag coefficient for a somewhat rough sphere*) g = 9.8;(*Gravity*) Manipulate[ Module[{sol, v0}, v0 = Sqrt[2 A/(ρm 4/3 π r^3 + B)]; sol = ParametricNDSolve[{ρm 4/ 3 π r^3 y''[t] == - ρm 4/3 π r^3 g - 1/2 ρa π r^2 Sqrt[x'[t]^2 + y'[t]^2] Cd y'[t], ρm 4/ 3 π r^3 x''[t] == -(1/2) ρa π r^2 Sqrt[ x'[t]^2 + y'[t]^2] Cd x'[t], x[0] == 0, y[0] == 0, x'[0] == v0 Cos[θ], y'[0] == v0 Sin[θ]}, {x[t], y[t]}, {t, 0, 100}, { r}]; FindMaximum[(((x[t] /. sol)[ R]) /. {FindRoot[(y[t] /. sol) [R], { t, 40, 0.05, 50}]})[[1]], {R, 0.05}] ], {{θ, π/4}, 0, π/2}, {{A, 149, "A [Throwing Ability] (J)"}, 0.0001, 200, Appearance -> {"Open"}}, {{B, 0.5, "B [~1/3 arm mass] (kg)"}, 0.0001, 2, Appearance -> {"Open"}}, Item[StringForm[ "\!$$\*SubscriptBox[\(ρ$$, $$m$$]\) = \ kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \ \!$$\*SubscriptBox[\(ρ$$, $$a$$]\) = \ kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \!$$\*SubscriptBox[\(C$$, \ $$d$$]\) = ", ρm, ρa // StandardForm, Cd // StandardForm]]]  • Parameter v0 isn't defined! Commented Oct 13, 2020 at 8:18 • Define a function for the root search outside of FindMaximum e.g.: funy[R_] := (((x[t] /. sol)[ R]) /. {FindRoot[{Evaluate[(y[t] /. sol)[R]] == 0}, {t, 40, 0.05, 50}]})[[1]]; Then you can print funy[R] and see what is does. Then you will see that the root search failes because it can not find a root in the interval given. Commented Oct 13, 2020 at 9:41 • Use WhenEvent  and Reap  and Sow[x[t]]  to stop integration of NDSolve at y[t]==0 and to get the x[t] at that t as the last value in the reap list. Plot the x values you then get as a function of r with the other parameters varied as you suggested with Manipulate. Commented Oct 13, 2020 at 10:36 • @Daniel I did, that works fine. Like I said it Just Works when I plot it Commented Oct 13, 2020 at 14:27 • Hi, if I do as you said, replacing FindMaximium by Plot and changing the range, I get the same erro message: "FindRoot::nlnum: The function value {ParametricFunction[1,InternalBag[<1>],2,1,{{r$6738},SystemUtilitiesHashTable[<2>],{},{},{1},{Automatic,0,0}},{NDSolvebase\$6747,NDSolveNDSolveParametricFunction[0,{ParametricNDSolve,InternalBag[<2>],None,ParametricNDSolve},<<6>>,{},All]}][R]} is not a list of numbers with dimensions {1} at {t} = {40.}. " Commented Oct 13, 2020 at 15:11

There seems to be two problem. FindMaximum has the attribute "HoldAll" and this prevents the replacements. Plot does not have this attribute, what eliminates the problem. What exactly is going on, I was not able to find out. But when you declare a function with the replacements outside of FindMaximum it works. The second problem is with the start value of R. It is too far off and the algorithm goes to the wrong side.

I post the changed code below.

ρm = 1600;(*Density of rock*)ρa = 1.225;(*Density of \
air*)Cd = 0.5;(*Typical drag coefficient for a somewhat rough \
sphere*)g = 9.8;(*Gravity*)Manipulate[
Module[{sol, v0}, v0 = Sqrt[2 A/(ρm 4/3 π r^3 + B)];
sol = ParametricNDSolve[{ρm 4/3 π r^3 y''[
t] == -ρm 4/3 π r^3 g -
1/2 ρa π r^2 Sqrt[x'[t]^2 + y'[t]^2] Cd y'[
t], ρm 4/3 π r^3 x''[
t] == -(1/2) ρa π r^2 Sqrt[x'[t]^2 + y'[t]^2] Cd x'[
t], x[0] == 0, y[0] == 0, x'[0] == v0 Cos[θ],
y'[0] == v0 Sin[θ]}, {x[t], y[t]}, {t, 0, 100}, {r}];
fun[R_ /;
NumericQ[
R]] := (((x[t] /. sol)[
R]) /. {FindRoot[(y[t] /. sol)[R], {t, 40, 0.05,
50}]})[[1]] /. {θ -> π/4, A -> 149, B -> 0.5};
FindMaximum[fun[R], {R, 0.01}]
], {{θ, π/4},
0, π/2}, {{A, 149, "A [Throwing Ability] (J)"}, 0.0001, 200,
Appearance -> {"Open"}}, {{B, 0.5, "B [~1/3 arm mass] (kg)"},
0.0001, 2, Appearance -> {"Open"}},
Item[StringForm[
"\!$$\*SubscriptBox[\(ρ$$, $$m$$]\) = \
kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \
\!$$\*SubscriptBox[\(ρ$$, $$a$$]\) = \
kg/\!$$\*SuperscriptBox[\(m$$, $$3$$]\), \!$$\*SubscriptBox[\(C$$, \
$$d$$]\) = ", ρm, ρa // StandardForm,
Cd // StandardForm]], TrackedSymbols -> {θ, A, B}]