import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


def solveit(n,y0):
    def exam(y, x):
        theta, omega = y
        dydx = [omega, - (2.0/x)*omega - theta**n]
        return dydx

    x = np.linspace(0.1, 10, 100)

    #call integrator
    sol = odeint(exam, y0, x)

    plt.plot(x, sol[:, 0], label='For n = %s,y0=(%s,%s)'%(n,y0[0],y0[1]))


ys= [[1.0, 0.0],[1.2, 0.2],[1.3, 0.3]]

fig = plt.figure()
for y_ in ys:
    solveit(1.,y_)

plt.legend(loc='best')
plt.grid()
plt.show()