Plotting a solution and its derivative, of a first order ODE

I have this code to solve a simple first order ODE using odeint. I managed to plot the solution y(r), but I also want to plot the derivative y'= dy/dr. I know y' it is given by f(y,r), but I'm not sure how to call the function with the output of the integration. Thank you in advance.

    from math import sqrt
    from numpy import zeros,linspace,array
    from scipy.integrate import odeint
    import matplotlib.pylab as plt

    def f(y,r):
        f = zeros(1)
        f[0] = -(y[0]*(y[0]-1.0)/r)-y[0]*(2/r+\
        ((r/m)/(1-r**2/m))*(2*sqrt(1-r**2/m)-k)/(k-sqrt(1-r**2/m)))\
        -(1/(1-r**2/m))*(-l*(l+1)/r+\
         (3*r/m)*(k+2*sqrt(1-r**2/m))/(k-sqrt(1-r**2/m)))\
        +((4*r**3)/((m**2)*(1-r**2/m)))*(1/(k-sqrt(1-r**2/m))**2)
        return f

    m = 1.15            
    k = 3*sqrt(1-1/m)
    l = 2.0
    r = 1.0e-10                         
    rf = 1.0                         

    rspan = linspace(r,rf,1000)
    y0 = array([l])
    sol = odeint(f,y0,rspan)
    plt.plot(rspan,sol,'k:',lw=1.5)

Solution

From odeint documentation:

For new code, use scipy.integrate.solve_ivp to solve a differential equation.

I have modified your code in the following manner and obtained the figure below.

import matplotlib.pyplot as plt
from numpy import gradient, squeeze, sqrt
from scipy.integrate import solve_ivp


def fun(t, y):
    l = 2
    m = 1.15
    k = 3 * sqrt(1 - 1 / m)
    return (-y * (y - 1) / t - y * (2 / t + t / m / (1 - t ** 2 / m)
                                    * (2 * sqrt(1 - t ** 2 / m) - k)
                                    / (k - sqrt(1 - t ** 2 / m)))
            - 1 / (1 - t ** 2 / m) * (-l * (l + 1) / t + 3 * t / m
                                      * (k + 2 * sqrt(1 - t ** 2 / m))
                                      / (k - sqrt(1 - t ** 2 / m)))
            + 4 * t ** 3 / m ** 2 / (1 - t ** 2 / m)
            / (k - sqrt(1 - t ** 2 / m)) ** 2)


sol = solve_ivp(fun, t_span=[1e-10, 1], y0=[2], method='BDF',
                dense_output=True)
if sol.success is True:
    print(sol.t.shape, sol.y.shape)
    plt.plot(sol.t, squeeze(sol.y), color='xkcd:avocado',
             label='Scipy Solution')
    plt.plot(sol.t, fun(sol.t, squeeze(sol.y)), linestyle='dashed',
             color='xkcd:purple', label='Derivative Using the Function')
    plt.plot(sol.t, gradient(squeeze(sol.y), sol.t), linestyle='dotted',
             color='xkcd:bright orange', label='Derivative Using Numpy')
    plt.legend()
    plt.tight_layout()
    plt.savefig('so.png', bbox_inches='tight', dpi=300)
    plt.show()

To address the comment about squeeze, please see below (extracted from scipy.integrate.solve_ivp):

where n is defined according to: