How to trace the curve of max points given sets of points

Question

I was working with Runge Kutta 4 and a quadratic friction. I was running some tests to see for which angle I would have the max range (it's Pi/4 as it would be for no friction). Anyway I got this image

First Picture

And I got that from this code (you don't have to read it, just here if anyone's interested):

#include <iostream>
#include <cmath>
#include <iomanip>
#include <fstream>

void RK4Step(double t, double *Y, void (*RHS)(double, double*, double*), double dt, int neq);
void RHS(double S, double *Y, double *R);

double B=4.e-5, g=9.81, m=1., h=1e-3, v=5;
int Neq=4;

int main() {

    std::cout<<std::setprecision(10);

    double th=0.;
    double Y[4];
   
    std::ofstream data; 

    data.open("Test.dat"); //Open file
    int Np=1e3;

    for (int i=0; i<17; i++) {
    
        Y[0]=0; Y[1]=0; Y[2]=v*cos(th); Y[3]=v*sin(th); //Reassign the initial values
    
        for (int j=0; j<Np; j++) { //For every angle keep trace of the projectile
        
            RK4Step(j*h, Y, RHS, h, Neq);
            data<<Y[0]<<"  "<<Y[1]<<std::endl;
        
        }
    
        data<<std::endl<<std::endl;
        th+=M_PI*0.03125; //Add to the angle
    }

    data.close(); //Close file.
}

void RK4Step(double t, double *Y, void (*RHS)(double, double*, double*), double dt, int neq) {

    double* Y1=new double[neq]; 
    double* k1=new double[neq];
    double* k2=new double[neq];
    double* k3=new double[neq];
    double* k4=new double[neq];

    RHS(t, Y, k1);

    for (int i = 0; i < neq; i++){
        Y1[i] = Y[i] + 0.5*dt*k1[i];
    }

    RHS(t + 0.5*dt, Y1, k2);

    for (int i = 0; i < neq; i++){
        Y1[i] = Y[i] + 0.5*dt*k2[i];
    }

    RHS(t + 0.5*dt, Y1, k3);

    for (int i = 0; i < neq; i++){
        Y1[i] = Y[i] + dt*k3[i];
    }

    RHS(t + dt, Y1, k4);

    for (int i = 0; i < neq; i++){
        Y[i] += dt*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i])*0.5/3.;
    }


    delete[] Y1; //Free up the memory
    delete[] k1;
    delete[] k2;
    delete[] k3;
    delete[] k4;

}

void RHS(double S, double *Y, double *R) {

     double vt=sqrt(Y[2]*Y[2]+Y[3]*Y[3]);

     R[0]=Y[2];
     R[2]=-B*Y[0]*vt;

     R[1]=Y[3];
     R[3]=-B*Y[1]*vt-g;

}

And plotted with gnuplot.

I was wondering is there any way to extrapolate this curve (the thick black one) from the data?

Second Picture

Maybe with an alghorithm that selects the maximun value of y for every point that has the around the same x, but definetely it would have to make big jumps for the first few poits (for example jumping from the peak of Pi/2 to the peak of 15Pi/32 without falling down) and then take points with a shorter intervalto have a smoother curve. But I wouldn't even know where to start. Or maybe I could directly do that with gnuplot. Any tips? Google searches didn't prove useful.