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

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()

• If your code looks like now after inserting the "code fences", then you have a simple indentation problem. Your main program code is, by indentation, still part of the function f, but without effect as after the unconditional return statement. Thus, fix the indentation and proceed to more interesting errors. Commented Oct 19, 2022 at 7:39
• @Lutz Lehmann, Thank you for your insightful comments. It took me quite a lot to understand indentation and finally make the correction only to finally see other "more interesting errors" as@ Lutz predicted. Also, the code is now showing like " Who is this 'p' and what is it doing here! as pointed out by Laurent, Besides, some new terms like 'operand' and 'scales' popped up. Learning to code is FUN indeed ! This time let me "try very hard" and yet still irritates me over and over I shall be coming back here shouting for help. Thank you Sirs, Commented Oct 19, 2022 at 17:53
• The component order of the input has to be the same in the output of f. The output is only the derivatives of the input components, nothing more. If you need the magnetic field components externally, use an extra function for the magnetic field, and call only this function to get the values of the field, in f and wherever else. Commented Oct 19, 2022 at 18:06
• Voila !!! Finally, the solutions are out. @LutzLehmann 'vielen Dank, mein Herr". With your help I have now learnt how to use different ODE solvers ( odeint, Rk4 and now solve-ivp. ). Please accept My arduous gratitude, respect and appreciation. Once I learnt enough I too shall always give time for enthusiatics learners. Commented Oct 19, 2022 at 18:44

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