from matplotlib import rcdefaults
rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'   
import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline   
sns.set()

from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook

import numpy as np

Runge-Kutta 4-stage 4th-order method

def rk4(f,t0,y0,h,m):

    b21 = 1/2
    b31 = 0
    b32 = 1/2
    b41 = 0
    b42 = 0
    b43 = 1
    b51 = 1/6
    b52 = 1/3
    b53 = 1/3
    b54 = 1/6
        
    t = t0
    y = np.array(y0)

    ya = np.empty((len(y0),m+1))
    ya[:,0] = y
    ta = np.linspace(t0,t0+m*h,m+1)
     
    for k in range(m):
        t = t0 + k*h
        K1 = f(t,y)  
        y1 = y + h*b21*K1
        K2 = f(t+h*b21,y1)
        y2 = y + h*(b31*K1 + b32*K2)
        K3 = f(t+h*(b31    + b32   ), y2)
        y3 = y + h*(b41*K1 + b42*K2 + b43*K3)
        K4 = f(t+h*(b31    + b32    + b43   ), y3)
        y4 = y + h*(b51*K1 + b52*K2 + b53*K3 + b54*K4)
        y = y4
        ya[:,k+1] = y
        
    return ta,ya

def f(x,Y):
    global p,q,phi
    y,yp = Y
    return np.array([ yp, phi(x) -p(x)*yp - q(x)*y ])

Plotting $y$ vs. $x$ for several initial slopes

def p(x):   return 5*x**2
def q(x):   return 30+50*np.sin(30*x)
def phi(x): return 50*np.cosh(x)
x0,x1 = 0,1
alpha = 1
beta = 10
y0 = alpha
m = 200
h = (x1-x0)/m
#for sigma in [2, 7, 5, 16, .25,-2, -10,-100]:
for sigma in np.linspace(0,-120,13):
    ta,Ya = rk4(f,x0,[y0,sigma],h,m)
    plt.plot(ta,Ya[0,:])

3D plot in $y$, $y^\prime$, $x$ space

fig = plt.figure(figsize=(10,7))
ax = fig.gca(projection='3d')
n=20
ya3 = np.empty((2,m+1,n))
for k,yp0 in enumerate( np.linspace(-30,30,n) ):
    ta,ya = rk4(f,0,[y0,yp0],1/m,m)
    ax.plot(ta,ya[0],ya[1],alpha = 0.95)
    ya3[:,:,k] = ya
ax.plot([0]*n,ya3[0, 0,:],ya3[1, 0,:],color='g')
ax.plot([1]*n,ya3[0,-1,:],ya3[1,-1,:],color='r')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('y \'')

Text(0.5, 0, "y '")