solve_ivp not giving out any output and no error while souple 3 coupled 2nd order ODES

Question

Please, someone tell me what is wrong in my code it does not give any outputs ( No plot nor print). The code is as below:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi

def f(t,p):
# assigning each ODE to a vector element
    r,θ,ϕ,x,y,z = p
    
    # constants
    Ω=9.74e-3
    B_θ=-8.6e-6*sin(θ)
    B_r=25893.2e-9*cos(θ)
    β=-9.36e-10
 
    # defining the ODEs
    dr_dt = x
    dx_dt = r*(y**2 + (z+Ω)**2 * sin(θ)**2 - β*z*sin(θ)*B_θ)
    dθ_dt = y
    dy_dt = (-2*x*y-r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r
    dϕ_dt = z
    dz_dt = (-2*x*(z+Ω)*sin(θ)-2*r*y*(z+Ω)*cos(θ)+β*(x*B_θ-r*y*B_r))/(r*sin(θ))

    return np.array([B_θ,B_r,dr_dt, dx_dt, dθ_dt, dy_dt, dϕ_dt, dz_dt])

    # time window
    t_span = (0, 100)
    t_eval = np.linspace(0,100, 200)
        
    # initial conditions
    p0 = np.array([0.7e+8,0.5*pi,0,0,0,0])
      
   # (5) Solve IVP
    sol = solve_ivp(f, t_span, p0, t_eval=t_eval)
    print (sol)
    plt.plot(sol.t, p[2])
    plt.show()

Answer 1

It seems you have not tried very hard to investigate the issue... I am pretty sure running your code raises an error which should guide you towards the solution.

You should plot something like sol.y[i,:], not p[2] which only exist within the scope of f.