Maximum distance traversed in an ideal projectile motion

Your method is not too bad at the start, except that you seem to have ignored (or at least not stated clearly) that $R$ depends on $theta$ as well, which is what makes this problem quite hard to solve. If I've understood you correctly, you would like to find the value of $theta$ (for a fixed speed $u$) that maximises the total length of the projectile in the air. In this case, taking derivatives like $frac{text{d}x}{text{d}theta}$ is not sensible. The variables that you want to maximise with respect to are $a,b,$ and $c$, since you will be integrating over $x$!

Despite how convinced I initially was that this problem should have a simple analytical result, it seems to not be the case! As far as I can see, to truly solve it, you need to use numerical methods. If anyone knows of a better way, I'd be very interested. Let me explain what I did.

I decided to make the following assumptions:

1. The total velocity (a constant) is 1. This is not a problem, I've just chosen units in which $u=1$, which is perfectly acceptable.

2. I will only be varying $u_y$, given the above constraint. The value of $u_x$ will be fixed by $sqrt{1 - u_y^2}$.

As you've pointed out (but formulated slightly differently) the total length covered by the projectile is:

$$L = int_{0}^{2u_y/g}sqrt{left(frac{text{d}y}{text{d}t}right)^2 + left(frac{text{d}x}{text{d}t}right)^2} text{d}t$$

(In this case, I have chosen to parameterise the curve by the time $t$, which I integrate from $t=0$ to $t=2 u_y/g$, which can be easily shown to be the total time of flight. You could do it your way too.)

Using the fact that

begin{equation} begin{aligned} y &= u_y t - frac{1}{2}g t^2,\ x &= u_x t, end{aligned} end{equation}

it's easy to show that

$$L = int_{0}^{2u_y/g}sqrt{(u_y-gt)^2 + u_x^2} ,,text{d}t = int_{0}^{2u_y/g}sqrt{u^2 - 2u_y g t + g^2 t^2} ,,text{d}t.$$

At times like this, it's useful to "adimensionalise" the equation, so that the limits don't depend on $u_y$. We can define a "dimensionless" time $$tau = frac{g}{2u_y}t,$$ so that the integral becomes:

$$L = frac{2}{g} int_{0}^{1}u_ysqrt{u^2 - 4 u_y^2 tau + 4 u_y^2 tau^2} ,,text{d}tau,$$

which is quite a nasty integral to solve by hand. Perhaps the people over at Math.SE would be able to do it justice? I decided to use Mathematica to solve it.

I first integrated the function numerically and plotted the integral for different values of $u_y$ as shown below, and was surprised to find that $L$ did have a maximum value (my initial thought was perhaps that it didn't) for $u_y$ somewhere between 0.82 and 0.84.

Given this, I got Mathematica to integrate the function and found that

$$L = frac{1}{4}left( 2 u u_y + (u_y^2 - u^2) lnleft({frac{u - u_y}{u+u_y}}right)right).$$

There is nothing stopping us from using units where $u=1$ and therefore $u_y in (0,1)$, and in these units

$$L= frac{1}{4}left( 2 u_y + (u_y^2 - 1) lnleft({frac{1 - u_y}{1+u_y}}right)right).$$

Next, I attempted to maximise this as a function of $u_y$ by taking the derivative and equating it to zero, which leads to:

$$2 + u_y lnleft({frac{1 - u_y}{1 + u_y}}right) = 0.$$

This is a transcendental equation and as such not easily solvable. But it's not too hard to solve it numerically to find that $L$ is maximised when $$u_y = 0.833557,$$

which lies in the range we expected. This corresponds to an angle of $$theta = arctan{frac{u_y}{sqrt{1 - u_y^2}}} = 0.985516 text{ rad} approx 56.466^circ.$$