FindMaximum not working with ParametricNDSolve

Question

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]]]

Answer 1

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}]