Why is my Phobos not orbiting around my Mars? It keeps making a straight line

I don't get how to fix it. I commented out some bits because I can't seem to find the initial velocity of Phobos. Here's my code:

import math
import matplotlib.pyplot as plt

G = 6.6743*(10**-11)
MMars = 0.64169 * (10**24)     #kg
MPhobos = 10.6 * (10**15)      #kg
#Vxphobos = 2138 
#Vxphobos = 7696700 
Vxphobos = 0
#Vxphobos = 3520
Vyphobos = 4100 
mars = [100,100]
phobos = [mars[1]-(6*(10**3)),0]
t = 0
phobosunit = [0,0]

def distt(distx,disty):
    totaldist = math.sqrt(distx**2+disty**2)
    return totaldist 

while t < 57552: #time in seconds that phobos takes to orbit mars  #27552
    tchange = 300
    t += tchange

    phobosr = distt((mars[0]-phobos[0]),(mars[1]-phobos[1]))
    phobosdist = (math.sqrt(((distt(phobos[0],mars[0]))**2)+((distt(phobos[1],mars[1]))**2)))
    for i in range(2):
        phobosunit[i] = -(mars[i]-phobos[i])/phobosdist
    plt.quiver(phobos[0],phobos[1],phobos[0]+phobosunit[0],phobos[1]+phobosunit[1])
    axphobos = (G*MMars/(phobosr**2))*phobosunit[0]
    ayphobos = (G*MMars/(phobosr**2))*phobosunit[1]
    #axphobos = (G*MMars/((distt(phobos[0],mars[0]))**2))*phobosunit[0]
    #ayphobos = (G*MMars/((distt(phobos[1],mars[1]))**2))*phobosunit[1]
    Vxphobos += axphobos*tchange
    Vyphobos += ayphobos*tchange
    phobos[0] += Vxphobos*tchange
    phobos[1] += Vyphobos*tchange
    plt.plot(phobos[0],phobos[1],'go') 
    print(phobosunit)
    force.append(G*MMars*MPhobos/(phobosr**2))

#plt.xticks([-3*e,-2*e,-1*e,0,1*e,2*e,3*e])
#plt.yticks([-3*e,-2*e,-1*e,0,1*e,2*e,3*e])

e = 10**8
plt.plot(mars[0],mars[1],'ro')
plt.xticks([-1*e,0,1*e,2*e,3*e])
plt.yticks([-1*e,0,1*e,2*e,3*e])
plt.show()

I used linearization to get from the force to the position. I feel like I made a really stupid mistake, but I have absolutely no idea on where to start to fix it. I've stared at this code for much too long now lmao. Any ideas?

It's making a straight line instead of a circle/oval.

Your equations of motion look alright to me. I think your constants are what's wrong.

If you're using SI units, then you should start with a speed of Vyphobos = 2.138e3 when it's at a height of ~9500km. Take a look at the following constants I defined, and their sources:

Phobos Mars

In addition, I used vector math with numpy to make the calculations easier. The vector form of gravity is a = G M R / r**3, where R is the vector between the bodies, and r is its norm.

import numpy as np
from matplotlib import pyplot as plt

G = 6.6743          # Exponent (-11) is multiplied later, in the variable E
MMars = 0.64169     # kg (exponent 24)

GM = G*MMars
E = 10**(24-11) # Exponent for G*Mmars

Vxphobos = 0.
Vyphobos = 2.138e3        # Vel is 2.138km/s

Rphobos = 11266.7         # Radius of phobos, m
Rmars = 3389.5e3          # Radius of Mars, m


mars = np.array([0., 0.])      # Mars at 0, 0
phobos = mars - [9500e3, 0]    # Phobos at a orbital radius of 9500km

Vphobos = np.array([Vxphobos, Vyphobos]) # Phobos velocity vector

tchange = 1                              # Smaller timestep
t = 0                                    #
tmax = 57552                             #
n_steps = (tmax - t) // tchange          # Number of timesteps
i = 0                                    # 

Now, define an array to hold the results:

phobos_loc_vel = np.zeros((n_steps+1, 4))

# Set the initial position
phobos_loc_vel[0, :2] = phobos
phobos_loc_vel[0, 2:] = Vphobos

Do the timestepping

while i < n_steps:
    i += 1
    t += tchange

    R = mars - phobos         # Vector to Mars from phobos
    r = np.linalg.norm(R)     # Distance between Mars and phobos

    
    aphobos = GM*R/r**3 * E      # Acceleration
    phobos += Vphobos*tchange    # Timestep position
    Vphobos += aphobos*tchange   # Timestep velocity

    # Store timestep in array
    phobos_loc_vel[i, :2] = phobos
    phobos_loc_vel[i, 2:] = Vphobos

And to plot:

fig, ax = plt.subplots()

# Plot every 300 timesteps
plt.plot(phobos_loc_vel[::300, 0], phobos_loc_vel[::300, 1], '.g')

# Draw Mars
c = plt.Circle(mars, Rmars, color='r')
ax.add_patch(c)

# Set plot limits and aspect ratio
ax.set_xlim([-15000e3, 15000e3])
ax.set_ylim([-15000e3, 15000e3])
ax.set_aspect('equal')

Which gives the following plot: