Trajectory with Thrust

We need to solve the simultaneous differential equations $m\ddot{x} \; = \; \frac{mg\dot{x}}{2\sqrt{\dot{x}^2+\dot{y}^2}} \hspace{2cm} m\ddot{y} \; = \; \frac{mg\dot{y}}{2\sqrt{\dot{x}^2+\dot{y}^2}} - mg$ with $x(0) = y(0) = 0$ and $\dot{x}(0) = V\cos\theta$ , $\dot{y}(0) = V\sin\theta$ for some $V > 0$ . By rescaling the time and distance variables we can simplify the equations to $\ddot{x} \; = \; \frac{\dot{x}}{2\sqrt{\dot{x}^2+\dot{y}^2}} \hspace{2cm} \ddot{y} \; = \; \frac{\dot{y}}{2\sqrt{\dot{x}^2+\dot{y}^2}} - 1$ with $x(0) = y(0) = 0$ and $\dot{x}(0) = 5\cos\theta$ , $\dot{y}(0) = 5\sin\theta$ . We then need to find the value $T_\theta>0$ for which $y(T_\theta) = 0$ , and seek to maximize $R(\theta) \,=\, x(T_\theta)$ . Here is a plot of $R(\theta)$ .

This is subject to numerical analysis. Solving these equations numerically using Mathematica, we obtain that $R$ is maximized when $\theta = 1.062339$ , which corresponds to $\boxed{60.87^\circ}$ .