Haz un programa en python que simule un sistema planetario de n cuerpos. Este preguntara por consola el numero de cuerpos, su masa, velocidad y posición inicial. Este ejecutara la simulación 10000 pasos donde cada paso es un segundo. Finalmente mostrara las trayectorias de los cuerpos en diferentes colores usando matplotlib.
¡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(fEnter velocity in x direction of body {i+1}: )), 
                           float(input(fEnter velocity in y direction of body {i+1}: ))]
          positions[i] = [float(input(fEnter x position of body {i+1}: )), 
                          float(input(fEnter 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.
Lo he probado con 2 cuerpos y funciona, pero cuando quiero probar con 3 o mas no funciona, ¿Puede explicarme el por que? y ayudarme a ejecutar la simulación con 3 o mas cuerpos?