Can we simulate with mathematical calculation using blender and python?
I won't pretend to know all required physics and the good initial conditions, but that can be done using the Python script.
The idea is to setup a function (a handler) which is fired by Blender each time frame (current time) is changing. This is done by &quot;bpy.app.handlers.frame_change_pre.append( theFunction )&quot; which mean &quot;theFunction&quot; is called before the change of the frame.
The result is the following:

The code:
import bpy
from mathutils import Vector

#get first object and arbitrary location
o1 = bpy.data.objects['o1']
o1.location = Vector( (-5,0,0 ) )

#same for second one
o2 = bpy.data.objects['o2']
o2.location = Vector( (5,0,0 ) )

#define constants: G and masses
G = 6.6740831e-11
m1 = 1e+11
m2 = 1e+05

#initial speeds
v1 = Vector( ( 0, 0.01, 0 ) )
v2 = Vector( ( 0, -0.5, 0 ) )

#frame rate for calculation
frame_rate = 1 #24 (normally 24 but 1 is faster)
seconds_per_frame = 1 / frame_rate

def two_bodies(scene):
    #import globals in the function scope
    global o1, o2, G, m1, m2, v1, v2
    global seconds_per_frame
    
    #calc current direction between the objects
    direction = o1.location - o2.location
    
    #calc the squared distance 
    d_squared = direction.length_squared
    #keep the direction of the strength
    direction.normalize()
    
    #calc new speed vectors
    v1 = v1 - (direction * (G * m2 / d_squared) * seconds_per_frame)
    v2 = v2 + (direction * (G * m1 / d_squared) * seconds_per_frame)
    
    #calc new locations
    o1.location += v1 * seconds_per_frame
    o2.location += v2 * seconds_per_frame

#get rid of previously set frame_change_pre handlers (if not the same handler may be fired n times)
bpy.app.handlers.frame_change_pre.clear()
#install &quot;two_bodies&quot; as current handler
bpy.app.handlers.frame_change_pre.append(two_bodies)

The blend file is below. To use it:

In the text part, use the &quot;run script&quot; button in order to setup initial conditions and the handler itself
In the 3D view part, use AltA to run the animation