v=10
theta=.7
solution= NDSolve[{
x''[t]== -Sqrt[x'[t]^2+z'[t]^2] x'[t],
z''[t]== -Sqrt[x'[t]^2+z'[t]^2] z'[t] - 1,
  x[0]== 0,
  z[0]== 0,
 x'[0]== v Cos[theta],
 z'[0]== v Sin[theta]},{x,z},{t,0,4}]

ParametricPlot[Evaluate[{x[t],z[t]} /. First[solution]],{t,0,3.2},PlotRange->All]




<<Cd.m
Plot[cd[m],{m,.1,10}]

R=.05
m=15
v=1000
rho=1.1
c=330

theta=???
Note...you're going to use somethig like:

cd[Sqrt[x'[t]^2+z'[t]^2] /c]

in your NDSolve...why?

