Thus it appears that at about Pi/4 or near .75 radians we have a maximum range. To determine the value of theta which maximizes this range precisely we determine where range'[theta]=0, i.e. where the derivative with respect to theta is 0.

Input := 


derrange[theta_] = range'[theta];

Input := 


sol = FindRoot[derrange[theta]== 0,{theta, Pi/8}]

Output =


{theta -> 0.785398}

Input := 


range[theta]/.sol[[1]]//N

Output =


2812.5