def ode_solve(y):
    r=sqrt(y[1]^2 + y[2]^2 + y[3]^2)
    n=1/(r^3+r^2+r)

    dy[1] = y[4]
    dy[2] = y[5]
    dy[3] = y[6]
    dy[4] = -y[1]*n/r
    dy[5] = -y[2]*n/r
    dy[6] = -y[3]*n/r

    return dy

def o_slv(y, t):
    x, v = y.reshape(2, -1)

    rsq = (x**2).sum()
    r   = np.sqrt(rsq)

    a   = -3*(x/r)/(r + rsq + rsq*r)  # use (x/r) to get the normal vector

    dy = np.hstack((v, a))

    return dy

import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import odeint as ODEint

velocities = np.linspace(0.7, 1.3, 7)
times      = np.linspace(0, 10, 1001)

answers = []
for v in velocities:

    Y0 = np.array([1, 0, 0, 0, v, 0])

    answer, info = ODEint(o_slv, Y0, times, full_output=True)
    answers.append(answer)

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    for a in answers:
        x, y, z = a.T[:3]
        plt.plot(x, y)
    plt.plot([0], [0], 'ok')
    plt.subplot(2, 1, 2)
    for i, name  in enumerate(('x', 'y', 'z')):
        plt.subplot(6, 1, i+1+3)
        for a in answers:
            plt.plot(times, a.T[i])
        plt.title(name, fontsize=16)
    plt.show()