assistant: ¡Claro! Aquí tienes un programa en Python que simula un sistema planetario.

```python
import numpy as np
import matplotlib.pyplot as plt

G = 6.67430e-11 # gravitational constant

def update_positions(positions, velocities, masses, time_step):
    n_bodies = len(masses)
    forces = np.zeros((n_bodies, 2))

    # calculate forces
    for i in range(n_bodies):
        for j in range(n_bodies):
            if i != j:
                r = positions[j] - positions[i]
                f = G * masses[i] * masses[j] / np.linalg.norm(r)**2
                forces[i] += f * r / np.linalg.norm(r)

    # update positions and velocities
    acceleration = forces / masses
    velocities += acceleration  * time_step
    positions  += velocities * time_step

    return positions, velocities

def simulate(n_bodies):
    if n_bodies == 2:
      masses     = 100*np.array([1, 1])
      velocities = np.array([[-0.05, 0.5], [0.05, 0.5]])
      positions  = np.array([[0.5, 0], [-0.5, 0]])

    else:
      # ask for input
      masses = np.zeros(n_bodies)
      velocities = np.zeros((n_bodies, 2))
      positions = np.zeros((n_bodies, 2))

      for i in range(n_bodies):
          masses[i] = float(input(f"\nEnter mass of body {i+1}: "))
          velocities[i] = [float(input(f"Enter velocity in x direction of body {i+1}: ")), 
                           float(input(f"Enter velocity in y direction of body {i+1}: "))]
          positions[i] = [float(input(f"Enter x position of body {i+1}: ")), 
                          float(input(f"Enter y position of body {i+1}: "))]

    # run simulation
    time_step = 1.0 # in seconds
    n_steps = 10000
    positions_history = np.zeros((n_steps+1, n_bodies, 2))
    positions_history[0] = positions.copy()
    for i in range(n_steps):
        positions, velocities = update_positions(positions, velocities, masses, time_step)
        positions_history[i+1] = positions.copy()

    # set up plotting
    colors = plt.cm.rainbow(np.linspace(0, 1, n_bodies))
    fig, ax = plt.subplots()
    for j in range(n_bodies):
        ax.plot(positions_history[:,j,0], positions_history[:,j,1], color=colors[j], linewidth=0.5)
        ax.plot(positions_history[-1,j,0], positions_history[-1,j,1], color=colors[j], marker='o', markersize=3)

    # show the plot
    plt.show()

if __name__ == '__main__':
    simulate(n_bodies = int(input("Enter number of bodies: ")))
```

El programa funciona de la siguiente manera:

1. Primero definimos la constante gravitatoria G.

2. La función *update_positions* actualiza las posiciones y velocidades de los cuerpos en función de sus masas y fuerzas gravitatorias.

3. La función *simulate* pide al usuario los datos de entrada (masas, velocidades y posiciones de cada cuerpo), configura el trazado, ejecuta la simulación durante 10000 pasos (cada paso es un segundo) y traza las trayectorias de los cuerpos.

4. En el bloque ```python if __name__ == '__main__' ```: pedimos al usuario el número de cuerpos a simular y llamamos a la función simulate.

Para ejecutar el programa, guarde el código como un archivo Python (por ejemplo, n_body_simulation.py) y ejecútelo desde la línea de comandos. Se le pedirá que introduzca el número de cuerpos y, a continuación, para cada cuerpo, tendrá que introducir su masa, velocidad en la dirección x, velocidad en la dirección y, posición x y posición y. Después de ejecutar la simulación, el programa mostrará un gráfico de las trayectorias de los cuerpos.