Switching coordinate charts mid-integration with solve_across_charts

edit

I am trying to utilize solve_across_charts to avoid coordinate singularities by switching from spherical to cartesian coordinates when integrating geodesics in the Kerr spacetime. However the code just keeps infinitely running. Following the documentation, my code is (9.3 version),

M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')
X.<t,r,th,ph>=M.chart(r"t r:(2.011,+infinity) th:(0,pi):\theta ph:(0,2*pi):periodic\phi") 
C.<T,x,y,z>=M.chart()
C.add_restrictions(sqrt(x**2+y**2+z**2) > 2.011) 
X_to_C=X.transition_map(C,[t,r*sin(th)*cos(ph),r*sin(th)*sin(ph), r*cos(th)])
C_to_X=C.transition_map(X,[T,sqrt(x**2+y**2+z**2),atan(sqrt(x**2+y**2)/z),atan(y/x)])
var('m,b,a', domain='real')
m=1
g = M.metric()

rho2 = r^2 + (a*cos(th))^2
Delta = r^2 -2*m*r + a^2
g[0,0] = -(1-2*m*r/rho2)
g[0,3] = -2*a*m*r*sin(th)^2/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2

def g00(m,r,th,a):return g[0,0](m,r,th,a)
def g03(m,r,th,a):return g[0,3](m,r,th,a)
def g11(m,r,th,a):return g[1,1](m,r,th,a)
def g22(m,r,th,a):return g[2,2](m,r,th,a)
def g33(m,r,th,a):return g[3,3](m,r,th,a)
def D(m,r,th,a): return (g03(m,r,th,a))^2-g00(m,r,th,a)*g33(m,r,th,a)

def initial_vector(r0, b,al,ph0=0, E=1, inward=False):
print(b)
t0,th0=0,pi/2
L = -b*E
vth0 = al/(r0*r0)
vt0=(E*g33(m,r0,th0,a)+L*g03(m,r0,th0,a))/(D(m,r0,th0,a))
vr0=sqrt((((1/D(m,r0,th0,a))*((L^2)*g00(m,r0,th0,a)+2*E*L*g03(m,r0,th0,a)+(E^2)*g33(m,r0,th0,a)))
          -((vth0)^2)*g22(m,r0,th0,a))/g11(m,r0,th0,a))
if inward:

    vr0 = - vr0

vph0 = -(1/D(m,r0,th0,a))*(E*g03(m,r0,th0,a)+L*(g00(m,r0,th0,a)))
p0 = M((t0, r0, th0, ph0), name='p_0')

return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')


r0 = 100 
th0=pi/2
s = var('s') 

ph0=0
v0 = initial_vector(r0, b=0,al=0, ph0=ph0, inward=True)
geod = M.integrated_geodesic(g, (s, 0, 150), v0, across_charts=True)

sol = geod.solve_across_charts(step=0.1, method='odeint', parameters_values={a:0.5},
                                               verbose='False')
interp = geod.interpolate()

I think the issue is on the declaration of the transition maps, however I've tried everything I could think of with no luck. Any help would be much appreciated.