Second order coupled ODE using ODEINT

I am new to solving coupled ODEs with python, I am wondering if my approach is correct, currently this code outputs a graph that looks nothing like the expected output. These are the equations I am trying to solve:

And here is the code I am using (for the functions f_gr, f_sc_phi and f_gTheta you can just put any constant value)

import Radial as rd
import ScatteringAzimuthal as sa
import PolarComponent as pc
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


#gamma for now set to 1
g_mm = 1

def f(u,t):
    #y1 = thetadot ::  y2 = phidot  ::  y3 = cdot
    rho, theta, y1, phi, y2, c, y3 = u
    p = [y1, (pc.f_gTheta(theta,524.1+rho)/(c*np.cos(phi))-(g_mm*y1)+(2*y1*y2*np.tan(phi))-(2*y3*y1/c)), 
         y2, ((sa.f_sc_phi(theta,524.1+rho/c))-(g_mm*y2)-(2*y3*y2/c)-(np.sin(phi)*np.cos(phi)*y2**2)), 
         y3, (rd.f_gr(theta,524.1+rho)-(g_mm*y3)+(c*y2**2)+(c*(y1**2)*(np.cos(phi)**2))), phi]
    return p

time = np.linspace(0,10,100)
z2 = odeint(f,[0.1,np.pi/2,0.1,np.pi/2,0.1,0.1,0.1], time)

rhoPl = z2[:,0]
thetaPl = z2[:,1]
phiPl = z2[:,3]

'''
plt.plot(rhoPl,time)
plt.plot(thetaPl,time)
plt.plot(phiPl,time)
plt.show()

'''

x = rhoPl*np.sin(thetaPl)*np.cos(phiPl)
y = rhoPl*np.sin(thetaPl)*np.sin(phiPl)
z = rhoPl*np.cos(thetaPl)

plt.plot(x,time)
plt.plot(y,time)
plt.plot(z,time)
plt.show()

when I change the time from 0.1 to 5 I get an error:

ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.

Any ideas on how to improve this code or if my approach is completely incorrect?

Code for Radial.py

import numpy as np
from scipy.special import spherical_jn
from scipy.special import spherical_yn
import sympy as sp
import matplotlib.pyplot as plt
R_r = 5.6*10**(-5) 
l = 720
n_w = 1.326
#k = 524.5/R_r
X_r = 524.5



# R is constant r is changing  
def f_gr(theta,x):
    f = ((sp.sin(theta))**(2*l-2))*(1+(sp.cos(theta))**2)
    b = (spherical_jn(l,n_w*x)*spherical_jn(l,n_w*x,True))+(spherical_yn(l,n_w*x)*spherical_yn(l,n_w*x,True))
    c = (spherical_jn(l,n_w*X_r)*spherical_jn(l,n_w*X_r,True))+(spherical_yn(l,n_w*X_r)*spherical_yn(l,n_w*X_r,True))
    n = b/c
    f = f*n
    return f

Code for ScatteringAzimuthal.py

from scipy.special import spherical_jn, spherical_yn
import numpy as np
import matplotlib.pyplot as plt

l = 720
n_w = 1.326
n_p = 1.572
X_r = 524.5
R_r = 5.6*10**(-5)
R_p = 7.5*10**(-7)
k = X_r/R_r


def f_sc_phi(theta,x):
    f = (2/3)*(n_w**2)*((X_r**3)/x)*((R_p**3)/(R_r**3))*(((n_p**2)-(n_w**2))/((n_p**2)+(2*(n_w**2))))
    g = np.sin(theta)**(2*l-3)
    numerator = (l*(1+np.sin(theta))- np.cos(2*theta))\
                *((spherical_jn(l,n_w*x)*spherical_jn(l,n_w*x))+(spherical_yn(l,n_w*x)*spherical_yn(l,n_w*x)))
    denominator = ((spherical_jn(l,n_w*X_r)*spherical_jn(l,n_w*X_r,True))\
                   +(spherical_yn(l,n_w*X_r)*spherical_yn(l,n_w*X_r,True)))
    m = numerator/denominator
    final = f*g*m
    return final

And Code for PolarComponent.py

import numpy as np
from scipy.special import spherical_yn, spherical_jn
import matplotlib.pyplot as plt

l = 720
n_w = 1.326
X_r = 524.5 #this value is implemented in the ode file

#define dimensionless polar component
#X_r is radius, x is variable 
def f_gTheta(theta,x):
    bessel1 = (spherical_jn(l,n_w*x)*spherical_jn(l,n_w*x)) + \
              (spherical_yn(l,n_w*x)*spherical_yn(l,n_w*x))
    bessel2 = ((spherical_yn(l,n_w*X_r)*spherical_yn(l,n_w*X_r,True)) + \
              (spherical_yn(l,n_w*X_r)*spherical_yn(l,n_w*X_r,True)))*n_w*x
    bessels = bessel1/bessel2
    rest = (np.sin(theta)**(2*l-3))*((l-1)*(1+(np.cos(theta)**2)) \
                                     -((np.sin(theta)**2)*np.cos(theta)))
    final = rest*bessels
    return final