Off[General::spell];
Off[General::spell1]; ClearAll["Global`*"]

\(\(\( (*3*) \)\(
\)\(ϵ = 10\^\(-6\); 
    R = 3000; v0 = 200; g = 9.81; k = 0.0052; θ0 = 0.2; 
    Rnew = \(v0\^2*Sin[2*θ0]\)\/g; Δ = 0.2; 
    Label[start]; 
    sol1 = DSolve[{\(\(x'\)'\)[t] == \(-k\)*\(x'\)[t], \(\(y'\)'\)[
              t] == \(-k\)*\(y'\)[t] - g, \(x'\)[0] == 
            v0*Cos[θ0], \(y'\)[0] == v0*Sin[θ0], 
          x[0] == 0, y[0] == 0}, {x[t], y[t]}, t]; 
    T0 = \(2  v0*Sin[θ0]\)\/g; 
    sol2 = FindRoot[\((y[t] /. sol1[\([1]\)])\) == 0, {t, T0}]; 
    T = t /. sol2[\([1]\)]; Rold = Rnew; 
    Rnew = Evaluate[\(x[t] /. sol1[\([1]\)]\) /. sol2[\([1]\)]]; 
    If[\((Rold - R)\)*\((R - Rnew)\) > 
        0, \ Δ = Δ/2]; 
    Print["Rold= ", Rold, ";  Rnew=", 
      Rnew, ";  θ0= ", θ0, ";  Δ= ", Δ, ";  T0= ", T0]; 
    If[Δ > ϵ, 
      If[Rnew > 
          R, \ θ0 = θ0 - Δ, θ0 = θ0 + Δ]; Goto[start]]; 
    ParametricPlot[Evaluate[{x[t], y[t]} /. sol1[\([1]\)]], {t, 0, T}, 
      PlotRange -> All]; Print["θ0 = ", θ0]\)\)\)

sol1 = 
      DSolve[{x''[t] == -k*x'[t], y''[
                t] == -k*y'[t] - g, x'[0] == 
              v0*Cos[θ0], y'[0] == v0*Sin[θ0], 
            x[0] == 0, y[0] == 0}, {x[t], y[t]}, t] /. {v0 -> 200, 
          g -> 9.81, k -> 0.0052}

sol2 = 
      DSolve[{x''[t] == -k*x'[t], y''[
                  t] == -k*y'[t] - g, x'[0] == 
                v0*Cos[θ0], y'[0] == v0*Sin[θ0], 
              x[0] == 0, y[0] == 0} /. {k -> 0}, {x[t], y[t]}, 
          t] /. {v0 -> 200, g -> 9.81, k -> 0}

sol3 = 
      Solve[{x[t] == 3000, y[t] == 0} /. 
            sol2[[1]] /. {t -> 
              Tt, θ0 -> θ0t}, {Tt, θ0t}]

θ0t = (θ0t /. sol3[[3]]); 
    Tt = (Tt /. sol3[[3]]); 
    sol4 = FindRoot[
        Evaluate[{x[t] == 3000, y[t] == 0} /. 
              sol1[[1]] /. {t -> T}], {T, 
          Tt}, {θ0, θ0t}]

sol5 = 
      NDSolve[{x''[t] == -0.0052 x'[t], y''[
                t] == -0.0052 y'[t] - 9.81, x'[0] == 
              200*Cos[θ0], y'[0] == 200*Sin[θ0], x[0] == 0, 
            y[0] == 0} /. {θ0 -> 0.448265}, {x[t], y[t]}, {t, 0, 
          17.4092}]; 
    plot1 = ParametricPlot[
        Evaluate[{x[t], y[t]} /. sol5[[1]]], {t, 0, 
          17.4092}]

sol6 = 
      NDSolve[{x''[t] == 0, y''[t] == -9.81, x'[0] == 
              200*Cos[θ0], y'[0] == 200*Sin[θ0], x[0] == 0, 
            y[0] == 0} /. {θ0 -> 0.448265}, {x[t], y[t]}, {t, 0, 
          17.4092}]; 
    plot2 = ParametricPlot[
        Evaluate[{x[t], y[t]} /. sol6[[1]]], {t, 0, 17.4092}, 
        PlotStyle  ->  { Dashing[{0.03,  .03}]}]

Show[plot1, plot2]