Simulando el sistema solar con la ley de Newton

Hice una simulación en C++ con la ley de Newton y la probé comparando las posiciones de los planetas con la posición de la calculadora del sistema solar Don Cross (que convertí de JavaScript a C++)
http://cosinekitty.com/solar_system.html

Lo que hago es cada paso de tiempo (generalmente 1 segundo, pero el paso de 0,2 segundos es muy similar al paso de 10 segundos):

1. Calcular aceleración ($=$newton se forza$\times$tiempo delta)
   
2. Velocidad de actualización y posiciones
   
3. Compare las posiciones con los resultados de los algoritmos de la calculadora solar Don Cross.
   Pero después de 10 días de simulación, obtengo esta desviación de distancia (a la calculadora de Don Cross):
   

$\mathrm{Mercury} \ 4498.7 \ \mathrm{km}$

$\mathrm{Venus \ X} \ 1939.8 \ \mathrm{km}$

$\mathrm{Earth \ X} \ 10614.6 \ \mathrm{km}$

$\mathrm{Moon \ X} \ 7800.2 \ \mathrm{km}$

$\mathrm{Mars \ X} \ 445.2 \ \mathrm{km}$

$\mathrm{Ceres \ X} \ 129.5 \ \mathrm{km}$

$\mathrm{Pallas \ X} \ 432.4 \ \mathrm{km}$

$\mathrm{Juno \ X} \ 151.4 \ \mathrm{km}$

$\mathrm{Vesta \ X} \ 157.6 \ \mathrm{km}$

$\mathrm{Ida \ X} \ 73.6 \ \mathrm{km}$

$\mathrm{Gaspra} \ 455.3 \ \mathrm{km}$

$\mathrm{9P/T1} \ 241.5 \ \mathrm{km}$

$\mathrm{19P/B} \ 402.7 \ \mathrm{km}$

$\mathrm{67P/C-G} \ 533.2 \ \mathrm{km}$

$\mathrm{81P/W2} \ 110.7 \ \mathrm{km}$

$\mathrm{Jupiter} \ 172.3 \ \mathrm{km}$

$\mathrm{Saturn} \ 261.2 \ \mathrm{km}$

$\mathrm{Uranus} \ 71.4 \ \mathrm{km}$

$\mathrm{Neptune} \ 31.3 \ \mathrm{km}$

$\mathrm{Pluto \ X \ } \ 45.7 \ \mathrm{km}$

Como puede ver, algunos planetas tienen pequeñas desviaciones y otros más grandes, por lo que mi pregunta es: ¿Puede ser exacto el de Newton? o la calculadora del sistema solar Don Cross no lo es? ¿O hay materia negra en esa región? ¿O qué más?

void CGravitator::CalcAceleration(double timeseconds){
unsigned int i,j,iend;
if (sunStatic)iend=m_np-1;
else iend=m_np;
for (i = 0; i < iend; i++) {
        m_planetas[i].aceleration.set(0,0,0);
        CVector3 totalGravitationalForce;                                       
        // Loop through all bodies in the universe to calculate their gravitational pull on the object (TODO: Ignore very far away objects for better performance)


     for (j = 0; j < m_np; j++) {
            if (i == j) continue; // No need to calculate the gravitational pull of a body on itself as it will always be 0.                                
            double distancia =CVector3::Distancia(m_planetas[i].pos,m_planetas[j].pos);
            double force = KGNEWTON * m_planetas[i].masa * m_planetas[j].masa / pow(distancia, 2);
            CVector3 forceDirection = CVector3::Normalize(m_planetas[j].pos - m_planetas[i].pos);
            
            totalGravitationalForce += forceDirection * force;
        }
            CVector3 incspeed = totalGravitationalForce / m_planetas[i].masa ;          
        m_planetas[i].aceleration += incspeed * timeseconds;
        
    
}