Why does my planetary orbit simulation produce a straight line instead of an elliptical orbit?

I'm trying to simulate the orbit of a planet using the compute_orbit method, but when I plot the resulting positions, I get a straight line instead of an expected elliptical orbit. Below are the relevant parts of my code.

Code Snippet

def get_initial_conditions(self, planet_name):
    planet_id = self.PLANETS[planet_name]
    obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
    eph = obj.vectors()
    position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
    velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   
    scale_factor_position = 1  
    scale_factor_velocity = 1  
    return {
        "position": np.array(position) / scale_factor_position,
        "velocity": np.array(velocity) / scale_factor_velocity
    }

def compute_orbit(self, central_mass=1.989e30, dt=10000, total_time=31536000):
    initial_conditions_data = self.get_initial_conditions(self.planet_name)
    initial_conditions = [
        initial_conditions_data['position'][0], initial_conditions_data['velocity'][0],
        initial_conditions_data['position'][1], initial_conditions_data['velocity'][1],
        initial_conditions_data['position'][2], initial_conditions_data['velocity'][2]
    ]

    def f(t, state):
        x, vx, y, vy, z, vz = state
        r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
        G = 6.67430e-11
        Fx = -G * central_mass * self.mass * x / r**3
        Fy = -G * central_mass * self.mass * y / r**3
        Fz = -G * central_mass * self.mass * z / r**3
        return [vx, Fx / self.mass, vy, Fy / self.mass, vz, Fz / self.mass]

    t_span = (0, total_time)
    t_eval = np.arange(0, total_time, dt)  
    sol = scipy.integrate.solve_ivp(
        f, 
        t_span, 
        initial_conditions, 
        t_eval=t_eval, 
        rtol=1e-3, 
        atol=1e-6
    )
    
    x = sol.y[0]  
    y = sol.y[2]  
    z = sol.y[4]  
    positions = np.column_stack((x, y, z))
    return positions

i have tried experimenting with different initial conditions but i always get straight lines plotting the data , below is a minimal reproduceable example

Code Snippet

from astroquery.jplhorizons import Horizons
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def get_initial_conditions(planet_id):
        obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
        eph = obj.vectors()
        position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
        velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   
        scale_factor_position = 1  
        scale_factor_velocity = 1  
        return {
            "position": np.array(position) / scale_factor_position,
            "velocity": np.array(velocity) / scale_factor_velocity
        }

def compute_orbit(central_mass=1.989e30,rotating_mass=5.972e24, dt=10000, total_time=31536000):
        initial_conditions_data = get_initial_conditions(399)
        initial_conditions = [
            initial_conditions_data['position'][0], initial_conditions_data['velocity'][0],
            initial_conditions_data['position'][1], initial_conditions_data['velocity'][1],
            initial_conditions_data['position'][2], initial_conditions_data['velocity'][2]
                            ]

        def f(t, state):
            x, vx, y, vy, z, vz = state
            r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
            G = 6.67430e-11
            Fx = -G * central_mass * rotating_mass * x / r**3
            Fy = -G * central_mass * rotating_mass * y / r**3
            Fz = -G * central_mass * rotating_mass * z / r**3
            return [vx, Fx / rotating_mass, vy, Fy / rotating_mass, vz, Fz / rotating_mass]


        t_span = (0, total_time)
        t_eval = np.arange(0, total_time, dt)  
        sol = scipy.integrate.solve_ivp(
        f, 
        t_span, 
        initial_conditions, 
        t_eval=t_eval, 
        rtol=1e-3, 
        atol=1e-6)
        x = sol.y[0]  
        y = sol.y[2]  
        z = sol.y[4]  
        positions = np.column_stack((x, y, z))
        return positions

def plot_orbit(positions):
    x = positions[:, 0]
    y = positions[:, 1]
    z = positions[:, 2]
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(x, y, z, label='Orbit Path', color='blue')
    ax.set_xlabel('X Position (m)')
    ax.set_ylabel('Y Position (m)')
    ax.set_zlabel('Z Position (m)')
    ax.set_title('Planet Orbit Simulation')
    ax.scatter(0, 0, 0, color='yellow', s=100, label='Central Mass (Sun)')
    ax.legend()
    ax.set_box_aspect([1, 1, 1]) 
    plt.show()


positions = compute_orbit()
plot_orbit(positions)

This problem has mismatched units.

In the following code:

    obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
    eph = obj.vectors()
    position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
    velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   

The obj.vectors() query does not return output in meters and meters per second. As the documentation describes, it returns output in AU and AU per day.

Column Name	Definition
...	...
x	x-component of position vector (float, au, X)
vx	x-component of velocity vector (float, au/d, VX)
...	...

Source.

Later on, you're using a gravitational constant which is only appropriate for units of kilogram, meter, and second.

G = 6.67430e-11

The simplest way to fix this is to convert the units that Horizons uses into meters and meters per second. Alternatively, you could change G to a value suitable for AU and days.

def get_initial_conditions(planet_id):
        obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
        eph = obj.vectors()
        meters_per_au = 1.496e+11
        seconds_per_day = 86400
        position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) * meters_per_au
        velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]]) * meters_per_au / seconds_per_day
        scale_factor_position = 1  
        scale_factor_velocity = 1  
        return {
            "position": np.array(position) / scale_factor_position,
            "velocity": np.array(velocity) / scale_factor_velocity
        }

You should get the following result.