Help with Blender animation type thing. (New to blender)

So I made a 2D physics simulator in python using pygame which simulated a solar system. I wanted to make a 3-D version using blender but I have gotten stuck on how to actually make it run. With pygame I simply trapped it in a while loop and let it play out and closed it when I finished. It didn’t take me long to figure that blender doesn’t like that and freezes when stuck in a while loop.
My code works for creating and moving the objects but I need a way to play out the animation. I was thinking something like making it so that every frame of the animation the last part of my code runs. (or I could simply manually create the solar system and have a code which only causes gravity) I am very new to blender as today is my first day and I only picked it up to make this program. I would love for anyone to point me in the right direction.
Here is my code: (I apologize as there is probably a way to nicely format this which I am not seeing.)

import bpy

def distance(coord1, coord2):
return ((coord1[0] - coord2[0])**2 + (coord1[1] - coord2[1])**2 + (coord1[2] - coord2[2])**2)**0.5

def add(coord1, coord2):
return [coord1[0] + coord2[0], coord1[1] + coord2[1], coord1[2] + coord2[2]]

def normalized(coord1, coord2): #Starting at 1 and pointing towards 2.
d = distance(coord1, coord2)
return [(coord2[0] - coord1[0])/d, (coord2[1] - coord1[1])/d, (coord2[2] - coord1[2])/d]

def multiply(coord1, scalar):
return [coord1[0]*scalar, coord1[1]*scalar, coord1[2]*scalar]

Time = 0.1
N = 2
pi = 3.14159265359
pos = []
mass = []
vel = []
rad = []
den = []





vel.append((0,-0.15,0)) #Note that having VCM != 0 the system will drift.
G = 10

for i in range(N):
bpy.context.scene.cursor_location = pos[i]
bpy.ops.mesh.primitive_uv_sphere_add()[i].dimensions = (rad[i], rad[i], rad[i])
while True: Worried about trapping things in for loops in blender.
for i in range(N):
for j in range(N):
if i != j:
Force = mass[j]G/distance(pos[i], pos[j])**2 #Technically Force/Mass
vel[i] = add(vel[i], multiply(normalized(pos[i], pos[j]), Force
for i in range(N):
Pos =[i].location[i].location = add(pos[i], Pos)

Oh, thanks python was giving me issues with vector math, and I have it working now. Thanks!