Here is the equation I obtained, representing the horizontal range of a projectile in air resistance as a function of its initial release angle-
$$x=\frac{V_0\cos(\theta)}{B}-\frac{V_0\cos(\theta)}{B}\,\exp\left(\frac{BV_0\sin(\theta)}{g}\right)\exp\left(W\left(\frac{BV_0\sin(\theta)\exp\left(\frac{BV_0\sin(\theta)}{-g}\right)}{-g}\right)\right), $$
where:
$x=$ horizontal distance travelled by projectile
$V_0=$ initial velocity of projectile
$g=$ acceleration due to gravity
$\theta=$ initial release angle
I have assumed that the force of air resistance acting on the projectile is equal to a constant k multiplied by the velocity of the projectile.
So, $F=kv.$
In my equation, $B= k/m,$ where $m$ is the mass of the projectile.
Also, $W(\cdot)$ is the $W$ Lambert function.
I have been trying to differentiate the expression with respect to $\theta$ in order find an expression for the value of $\theta$ that would result in maximum horizontal range, assuming all other terms in the equation are constant.
How can this be done?
Starting from $\mathbf{F}=-k\,\mathbf{v}-mg\,\hat{\mathbf{j}}=m\dot{\mathbf{v}},$ and initial conditions $\mathbf{v}_0=v_0\langle\cos(\theta),\sin(\theta)\rangle$ and $\mathbf{x}_0=\langle 0,0\rangle,$ we obtain \begin{align*} m\ddot{x}&=-k\dot{x} \\ \ddot{x}&=-B\dot{x} \\ \dot{x}&=Ce^{-Bt} \\ \dot{x}&=v_0\cos(\theta)\,e^{-Bt} \\ x&=-\frac{v_0\cos(\theta)}{B}\,e^{-Bt}+C\\ x(t)&=\frac{v_0\cos(\theta)}{B}\,\left(1-e^{-Bt}\right),\;\text{and} \\ m\ddot{y}&=-k\dot{y}-mg \\ \ddot{y}&=-B\dot{y}-g \\ \int\frac{d \dot{y}}{B\dot{y}+g}&=-\int dt\\ \frac1B\int\frac{d \dot{y}}{\dot{y}+g/B}&=-t+C\\ \frac1B\,\ln|\dot{y}+g/B|&=-t+C\\ \ln|\dot{y}+g/B|&=-Bt+C, \quad\text{absorbing into constant}\\ \dot{y}+g/B&=Ce^{-Bt},\quad\text{switching to multiplicative constant}\\ v_0\sin(\theta)+\frac{g}{B}&=C\\ \dot{y}&=\left(v_0\sin(\theta)+\frac{g}{B}\right)e^{-Bt}-\frac{g}{B}\\ y&=-\left(\frac{v_0\sin(\theta)+\frac{g}{B}}{B}\right)e^{-Bt}-\frac{gt}{B}+C\\ y(t)&=\frac{v_0\sin(\theta)+mg/k}{B}\left(1-e^{-Bt}\right)-\frac{mgt}{k}. \end{align*} Solving the second equation $y(t)=0$ for $t$ yields $$t_{\text{range}}=\frac{mg W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+mg+k v_0 \sin (\theta )}{B mg}. $$ Plugging this back into $x$ yields $$x_{\text{range}}(\theta)=\frac{v_0 \cos (\theta )}{B} \left(1-\exp \left(- W\left(-e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (1+(k v_0/(mg)) \sin (\theta ))\right)+1+(k v_0/(mg)) \sin (\theta )\right)\right). $$ Now it is a fact that $$\dot{W}(t)=\frac{W(t)}{t(1+W(t))}. $$ Using this, we can differentiate this expression (it's quite a monster) to obtain $$x_{\text{range}}'(\theta)=-\frac{v_0 \left(mg W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+mg+k v_0 \sin (\theta )\right) \left((mg \sin (\theta )+k v_0) W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+\sin (\theta ) (mg+k v_0 \sin (\theta ))\right)}{B (mg+k v_0 \sin (\theta ))^2 \left(W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+1\right)}. $$ Solving $x_{\text{range}}'(\theta)=0$ for $\theta$ will be tricky. We can eliminate the denominator and just solve $$-v_0 \left(mg W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+mg+k v_0 \sin (\theta )\right) \left((mg \sin (\theta )+k v_0) W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+\sin (\theta ) (mg+k v_0 \sin (\theta ))\right) =0.$$ Mathematica does not immediately find a solution, but we note that this expression is factored. That is, we can solve either of \begin{align*} mg W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+mg+k v_0 \sin (\theta )&=0,\quad\text{or}\\ (mg \sin (\theta )+k v_0) W\left(-\frac{e^{-\frac{k v_0 \sin (\theta )}{mg}-1} (mg+k v_0 \sin (\theta ))}{mg}\right)+\sin (\theta ) (mg+k v_0 \sin (\theta ))&=0. \end{align*} The first one initially looks a bit more promising. Let $s=\sin(\theta),$ so that we have \begin{align*} mg W\left(-\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}\right)+mg+k v_0 s&=0 \\ W\left(-\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}\right)+1+(k v_0/(mg)) s&=0 \\ W\left(-\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}\right)&=-1-\frac{kv_0s}{mg} \\ -\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}&=-\left(\frac{mg+kv_0s}{mg}\right)\exp\left(-1-\frac{kv_0s}{mg}\right), \end{align*} which you can see is a tautology! That is, it conveys no information. Therefore, any solution must be from the other factor. Let us work on that (again, $s=\sin(\theta)$) \begin{align*} (mgs+k v_0) W\left(-\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}\right)+s (mg+k v_0 s)&=0 \\ (mg s+k v_0) W\left(-\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}\right)&=-s (mg+k v_0 s) \\ W\left(-\frac{e^{-\frac{k v_0 s}{mg}-1} (mg+k v_0 s)}{mg}\right)&=-s\,\frac{mg+kv_0s}{mgs+kv_0} \\ -\frac{e^{-\frac{k v_0 s-mg}{mg}} (mg+k v_0 s)}{mg}&=-s\,\frac{mg+kv_0s}{mgs+kv_0}\,\exp\left(-s\,\frac{mg+kv_0s}{mgs+kv_0}\right) \\ \exp\left(-\frac{k v_0 s-mg}{mg}\right) &=\frac{mgs}{mgs+kv_0}\,\exp\left(-s\,\frac{mg+kv_0s}{mgs+kv_0}\right). \end{align*} In canceling $mg+kv_0s,$ we are eliminating the solution $s=-mg/(kv_0),$ corresponding to $$\theta=\arcsin\left(-\frac{mg}{kv_0}\right). $$ A negative answer seems physically impossible, so we rule it out. We continue: \begin{align*} \exp\left(\frac{mg-k v_0 s}{mg}\right) &=\frac{mgs}{mgs+kv_0}\,\exp\left(-s\,\frac{mg+kv_0s}{mgs+kv_0}\right)\\ \frac{mg-k v_0 s}{mg}&=\ln\left(\frac{mgs}{mgs+kv_0}\right)-s\,\frac{mg+kv_0s}{mgs+kv_0}\\ \frac{g-B v_0 s}{g}&=\ln\left(\frac{gs}{gs+Bv_0}\right)-s\,\frac{g+Bv_0s}{gs+Bv_0}. \end{align*} Unfortunately, this is pretty much as far as we can go. $s$ appears in far too many places in this equation to think that there is an analytical solution. My recommendation is to take this equation, plug in the values you need for $m, g, k,$ and $v_0,$ and then solve it numerically.