Robust numerical computation of the 3D scalar potential field of the cubic Galileon gravity model at solar system scales