Ciao a tutti scrivo qui perchè avrei bisogno di un aiuto nel risolvere un problema che riguarda il metodo d'integrazione per le equazioni di un moto bidimensionale di un grave sottoposto ad un attrito quadratico (cioè quello dell'aria)

Questo è il codice, il cui unico scopo è capire se l'algoritmo d'integrazione funziona.

Il problema è che, prendendo come esempio un problema del libro, il risultato della gittata massima dovrebbe essere circa 59 m , mentre il programma mi da 58,14 metri
Ora, a me sembra un errore un po troppo consistente per dire che prendendo in considerazione un altro metodo d'integrazione possa darmi risultati migliori.

Sapreste darmi un parere su dove potrebbe celarsi il problema o se l'algoritmo è esatto e quindi non c'è nessun problema?

codice:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define g 9.81
#define PI 3.14159
#define A 0.0081666

struct variables {
  double x,y,v,vx,vy,alpha;
};

double AccelX(struct variables * p_XV);
double AccelY(struct variables * p_XV);

void rungeKutta2(struct variables * p_XV, double dt);

/*******************************************************************************************************************/

int main(int argc, char *argv[]) {
  int  i, numberOfSteps;
  double totalTime, dt;
  struct variables XV;
  FILE *fp;

  if (argc != 7) {
    fprintf(stderr,
	    "\nusage:   <x0> <y0> <v0> <alpha> <totalTime> <dt>\n\n");
    exit(1);
  }
  
  system("clear");
  
  XV.x = atof(argv[1]);
  XV.y = atof(argv[2]);
  XV.v= atof(argv[3]);
  XV.alpha=0.5*PI*atof(argv[4])/90.;
  totalTime = atof(argv[5]);
  dt = atof(argv[6]);
  XV.vx=XV.v*cos(XV.alpha);
  XV.vy=XV.v*sin(XV.alpha);
  
//  XV.vx=19.3;
//  XV.vy=23.0;
  
  printf("alpha= %lf   vx= %lf    vy= %lf\n\n",XV.alpha,XV.vx,XV.vy);
  getchar();
  numberOfSteps = (int)(totalTime / dt + 0.5);
    
  fp=fopen("RungeKutta.dat","w");
  fprintf(fp,"# Usiamo l'algoritmo di RungeKutta\n");
  fprintf(fp,"# 1:t  2:x  3:y  4:vx  5:vy\n");
  for (i = 0; i <= numberOfSteps; i++) {
     fprintf(fp,"%g %g %g %g %g\n", (double)i*dt, XV.x, XV.y,XV.vx,XV.vy);
     rungeKutta2(&XV,dt);
     if(XV.y<=0) i=numberOfSteps;
  }
  printf("t= %g   x= %g    y= %g    vx= %g     vy=  %g\n\n", (double)i*dt, XV.x, XV.y,XV.vx,XV.vy);
  
  fclose(fp);
  return EXIT_SUCCESS;
}

/*******************************************************************************************************************/

double AccelX(struct variables * p_XV) {
	return -A*p_XV->vx*sqrt(p_XV->vx*p_XV->vx+p_XV->vy*p_XV->vy);
}

double AccelY(struct variables * p_XV) {
     return -g-A*p_XV->vy*sqrt(p_XV->vx*p_XV->vx+p_XV->vy*p_XV->vy);
     
}
/*******************************************************************************************************************/

void rungeKutta2(struct variables * p_XV, double dt) {
  struct variables tmp;
  
  tmp.x = p_XV->x + 0.5 * p_XV->vx * dt;
  tmp.y = p_XV->y + 0.5 * p_XV->vy * dt;
  
  tmp.vx = p_XV->vx + 0.5 * AccelX(p_XV) * dt;
  tmp.vy = p_XV->vy + 0.5 * AccelY(p_XV) * dt;
  
  p_XV->x += tmp.vx * dt;
  p_XV->y += tmp.vy * dt;
  
  p_XV->vx += AccelX(&tmp) * dt;
  p_XV->vy += AccelY(&tmp) * dt;
}