Построение нескольких траекторий

Имеется проблема, необходимо построить модель движения первых пяти планет в солнечной системе, движение вокруг солнца, путем решения системы дифференциальных уравнений методом Рунге-Кутта 4 порядка. Результат работы модели - траектория.

Привожу код на примере вращения земли вокруг солнца:

import matplotlib.pyplot as plt
from numpy import zeros, linspace
from matplotlib.pyplot import style, figure, axes


# Функция f подготавливает массив, содержащий элементы вектор-функции,
# определяющей правую часть решаемой системы ОДУ
def f(m_sun, u, G):
    f = zeros(4)
    f[0] = u[2]
    f[1] = u[3]
    f[2] = - G * m_sun * u[0] / (u[0] ** 2 + u[1] ** 2) ** (3 / 2)
    f[3] = - G * m_sun * u[1] / (u[0] ** 2 + u[1] ** 2) ** (3 / 2)
    return f


# Определение входных данных задачи
# Определение входных данных задачи
t_0 = 0.
T = 365.25 * 24 * 60 * 60
x1_0 = 147098291 * 10 ** 3   # удаленность от солнца
x2_0 = 347098291 * 10 ** 3
x3_0 = 547098291 * 10 ** 3
y1_0 = 0.
y2_0 = 0.
y3_0 = 0.
v_x1_0 = 0.
v_x2_0 = 0.
v_x3_0 = 0.
v_y1_0 = 30.4 * 10 ** 3
v_y2_0 = 46.4 * 10 ** 3
v_y3_0 = 98.4 * 10 ** 3  # скорость движения по орбите
G = 6.674301515151515 * 10 ** (-11)
m_sun = 1.98847 * 10 ** 30

x_0 = [x1_0, x2_0, x3_0]
y_0 = [y1_0, y2_0, y3_0]
v_x_0 = [v_x1_0, v_x2_0, v_x3_0]
v_y_0 = [v_y1_0, v_y2_0, v_y3_0]

# Определение числа интервалов сетки,
# на которой будет искаться приближённое решение
M = 365

s = 4  # ERK4 - Рунге-Кутта 4 порядка для решения системы дифференциальных уравнений
# определяем требуемые коэффициенты для классической схемы 4 порядка
b = zeros(s); a = zeros((s,s)); c = zeros(s)
b[0] = 1/6; b[1] = 1/3; b[2] = 1/3; b[3] = 1/6
a[1,0] = 1/2; a[2,0] = 0.; a[2,1] = 1/2; a[3,0] = 0.; a[3,1] = 0.; a[3,2] = 1.
c[0] = 0.; c[1] = 1/2; c[2] = 1/2; c[3] = 1.0

# Определение сетки
dt = (T - t_0)/M
t = linspace(t_0,T,M + 1)

# Выделение памяти под массив сеточных значений решения системы ОДУ
# В строке с номером m этого массива хранятся сеточные значения решения,
# соответствующие моменту времени t_m
u = zeros((M + 1,4))

# Задание начальных условий
# (записываются строку с номером 0 массива u)
u[0] = [x_0, y_0, v_x_0, v_y_0]

# Реализация схемы ERKs
# (отметим, что во вспомогательных массивах κ, y_res и w
# размерность 4 соответствует числу компонент вектор-функции f)
for m in range(M):
    w = zeros((s,4))
    for p in range(s):
        k = zeros(4)
        for l in range(p):
            k = k + a[p, l]*w[l]
        w[p] = f(m_sun, u[m] + dt * k, G)
    y_res = zeros(4)
    for p in range(s):
        y_res = y_res + b[p]*w[p]
    u[m + 1] = u[m] + dt*y_res

# Отрисовка решения
fig = figure()
ax = axes(xlim=(-2*10**11,2*10**11), ylim=(-2*10**11,2*10**11))
ax.set_aspect('equal'); ax.set_xlabel('x'); ax.set_ylabel('y');
ax.plot(0,0,'yo',markersize=25)
ax.plot(u[:,0],u[:,1],'-g',markersize=5)
ax.plot(u[M,0],u[M,1], color='g', marker='o', markersize=7)
ax.set_title('Траектория движения Земли')
plt.show()

Вопрос заключается в следующем: каким образом, задать вектор начальных условий, для пяти планет, ибо сама система дифференциальных уравнений везед будет одинакова, чтобы получить 5 траекторий,т.е. для пяти планет? Т.е. я хочу сделать так, чтобы мог записать вектор начальных условйи для пяти планет и потом соответственно получить 5 траекторий...