```
%display latex
l=2
m=3
M = Manifold(4, 'M', structure='Lorentzian')
X.<t,r,th,ph> = M.chart('t r th ph')
U=(2*m*t+l^2-t^2)/(t^2+l^2)
g = M.metric()
g[0,0] = -1/U
g[1,1] = 4*l^2*U
g[1,3] = 4*l^2*U*cos(th)
g[2,2] = t^2+l^2
g[3,3] = 4*l^2*U*(cos(th))^2+(t^2+l^2)*(sin(th))^2
g.display()
R3 = Manifold(3, 'R^3', latex_name=r'\mathbb{R}^3')
X3.<x,y,z> = R3.chart()
to_R3 = M.diff_map(R3, {(X, X3): [r*sin(th)*cos(ph),
r*sin(th)*sin(ph), r*cos(th)]})
to_R3.display()
p0 = M.point((7, 30, 0, 0), name='p_0')
v0 = M.tangent_space(p0)((1, 53/12, 0, 0), name='v_0')
v0.display()
s = var('s')
geod = M.integrated_geodesic(g, (s, 0, 10), v0); geod
sol = geod.solve()
interp = geod.interpolate()
graph = geod.plot_integrated(chart=X3, mapping=to_R3, label_axes=True, plot_points=10)
graph += p0.plot(chart=X3, mapping=to_R3)
graph += sphere(color='grey')
show(graph)
lsoda-- at t (=r1), too much accuracy requested
for precision of machine.. see tolsf (=r2)
in above, r1 = 0.1000000000000D+00 r2 = NaN
/home/sc_serv/sage/local/var/lib/sage/venv-python3.10/lib/python3.10/site-packages/scipy/integrate/_odepack_py.py:247: ODEintWarning: Excess accuracy requested (tolerances too small). Run with full_output = 1 to get quantitative information.
warnings.warn(warning_msg, ODEintWarning)
lsoda-- at t (=r1), too much accuracy requested
for precision of machine.. see tolsf (=r2)
in above, r1 = 0.1000000000000D+00 r2 =
```

Using

one can check that sin(th) appears in the denominator of the last geodesic equation, so th=0 in the definition of p0, v0 is problematic.