In[]:=
sol=First@With[{g=9.81,b=0.5,theta=45Degree,v=20},NDSolve[{x''[t]==-bx'[t],y''[t]==-b*y'[t]-g,x[0]==0,y[0]==0,x'[0]==vCos[theta],y'[0]==vSin[theta]},{x,y},{t,0,10}]];​​
In[]:=
ParametricPlot[{x[t]/.sol,y[t]/.sol},{t,0,3}]
Out[]=
In[]:=
tfinal[theta_]:=First@With[{g=9.81,b=0.5,v=20},NDSolve[{x''[t]==-bx'[t],y''[t]==-by'[t]-g,x[0]==0,y[0]==0,x'[0]==vCos[thetaDegree],y'[0]==vSin[thetaDegree]},{x,y},{t,0,10}]];
In[]:=
ParametricPlot[Evaluate[Table[{x[t]/.tfinal[m],y[t]/.tfinal[m]},{m,20,80,10}]],{t,0,4},PlotRange->{Automatic,{0,15}},ImageSizeLarge]
Out[]=
flightTime=FindRoot[Evaluate[(y[t]/.tfinal[45])]==0,{t,3}](*Theta=45Degree*)
Out[]=
{t2.41038}
In[]:=
maxRange=Evaluate[x[t]/.tfinal[45]]/.flightTime(*Theta=45Degree*)
Out[]=
19.8093
In[]:=
data=Table[{theta,Evaluate[x[t]/.tfinal[theta]]/.FindRoot[Evaluate[(y[t]/.tfinal[theta])]==0,{t,2}]},{theta,20,60,1/10}];
In[]:=
ifun=Interpolation[data];Plot[ifun[theta],{theta,20,60}]
Out[]=
FindMaximum[ifun[theta],{theta,45}](*Maximumrangeis20.763795443462303[meters]forOnlytheta=35.46015884691499[degree]andfor:g=9.81,b=0.5,v=20*)
Out[]=
{20.7638,{theta35.4602}}