circular orbit weird with physics

Well, nothing you have experienced was unexpected or is a bug or anything like that. When you apply Newtons law of gravitation every frame you essentially performing dead reckoning in discrete time steps. This will never be "accurate" or "stable" over time since even the slightest error will accumulate.

Keep in mind that the real world does not "run" in discrete time intervals and also doesn't have any numerical precision limit on any involved figures. Also keep in mind that orbits are almost always elliptical and not circular.

If you want a stable orbit simulation you can not just use newton's law and integrate it over time with singe precision floating point values. You would need to specify your exact orbit and use absolute figures (like orbital elements)to specify your position on that orbit.

Note that almost all time discrete physics systems will have this issue since any change between two time steps is considered linear. Any non linear process will accumulate errors over time. This error gets smaller if you do smaller time steps, however the realworld does actual continuous integration. So in theory infinitely small timesteps. Though even when you have unlimited time to calculate infinite steps you would still accumulate errors due to number precision limits of floating point numbers.

ps: 2 units off when having a 2000 unit radius is only an error of 0.1% after how many iterations?

Which integration method are you using? I find that when using Euler integration, strange orbit behaviour occurs, like non-constant speed when in a circular orbit.

Here is some C# code for doing symplectic or Runge-Kutta integration:

 void proceed_euler(ref Vector3 pos, ref Vector3 vel, Vector3 ang_vel, float dt)
 {
     vel = vel + acceleration(pos, vel, ang_vel) * dt;
     pos = pos + vel * dt;
 }
 
 void proceed_symplectic4(ref Vector3 pos, ref Vector3 vel, Vector3 ang_vel, float dt)
 {
     float cr2 = Mathf.Pow(2.0f, 1.0f / 3.0f);
     
     float[] c = new float[4]
     {
         1.0f / (2.0f*(2.0f - cr2)),
         (1.0f - cr2) / (2.0f*(2.0f - cr2)),
         (1.0f - cr2) / (2.0f*(2.0f - cr2)),
         1.0f / (2.0f*(2.0f - cr2))
     };
     
     float[] d = new float[4]
     {
         1.0f / (2.0f - cr2),
         -cr2 / (2.0f - cr2),
         1.0f / (2.0f - cr2),
         0.0f
     };
     
     pos += vel * c[0] * dt;
     vel += acceleration(pos, vel, ang_vel) * d[0] * dt;
     
     pos += vel * c[1] * dt;
     vel += acceleration(pos, vel, ang_vel) * d[1] * dt;
     
     pos += vel * c[2] * dt;
     vel += acceleration(pos, vel, ang_vel) * d[2] * dt;
     
     pos += vel * c[3] * dt;
     //    vel += acceleration(pos, vel, ang_vel) * d[3] * dt; // last element d[3] is always 0
 }
 
 void proceed_rk4(ref Vector3 pos, ref Vector3 vel, Vector3 ang_vel, float dt)
 {
     const float one_sixth = 1.0f / 6.0f;
     
     Vector3 k1_velocity = vel;
     Vector3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
     Vector3 k2_velocity = vel + k1_acceleration * dt * 0.5f;
     Vector3 k2_acceleration = acceleration(pos + k1_velocity * dt * 0.5f, k2_velocity, ang_vel);
     Vector3 k3_velocity = vel + k2_acceleration * dt * 0.5f;
     Vector3 k3_acceleration = acceleration(pos + k2_velocity * dt * 0.5f, k3_velocity, ang_vel);
     Vector3 k4_velocity = vel + k3_acceleration * dt;
     Vector3 k4_acceleration = acceleration(pos + k3_velocity * dt, k4_velocity, ang_vel);
     
     vel += (k1_acceleration + (k2_acceleration + k3_acceleration) * 2.0f + k4_acceleration) * one_sixth * dt;
     pos += (k1_velocity + (k2_velocity + k3_velocity) * 2.0f + k4_velocity) * one_sixth * dt;
 }