StiffSystem or Singularity - a system of second order ODEs in the problem of geodesics

I would be extremely grateful for any help regarding the following code I wrote and the errors it produces. In this code I am investigating the behaviour of a massive particle trapped in the vicinity of the gravitational Bessel wave vertex.

I am really sorry for the lengthy expressions - I have no idea how to shorten it in any way that would help significantly. The table defined below is the right hand side of the geodesic equation. I have checked its correctness with my advisor, obtaining the same result as he did. The parameter "a" refers to the amplitude of the wave, while "b" and "c" are responsible for the wavevector quantities.

Remove["Global*"]
coords = {t, x, y, z}
n = Length[coords];

geodesic = {BesselJ[3,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] (-((a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) (x')^2)/(x^2 +
y^2)^(3/2)) + (2 a b^4 Sqrt[
b^2 + c^2] y (-3 x^2 + y^2) x' y')/(x^2 + y^2)^(3/
2) + (a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) (y')^2)/(x^2 +
y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b^4 Sqrt[b^2 + c^2] y (-3 x^2 + y^2) (x')^2)/(x^2 +
y^2)^(3/2) + (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(3/
2) - (a b^4 Sqrt[
b^2 + c^2] (-3 x^2 y + y^3) (y')^2)/(x^2 + y^2)^(3/2))) +
BesselJ[0,
b Sqrt[x^2 + y^2]] (4 a b Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^3 Cos[
Sqrt[b^2 + c^2] t - c z] x' z' +
4 a b Sqrt[b^2 + c^2] (c + Sqrt[b^2 + c^2])^3 Sin[
Sqrt[b^2 + c^2] t - c z] y' z') +
BesselJ[2,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (x^2 - y^2) x' z')/(x^2 +
y^2) + (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y y' z')/(x^2 + y^2)) -
Sin[
Sqrt[b^2 + c^2] t -
c z] (-((8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y x' z')/(x^2 +
y^2)) - (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (-x^2 +
y^2) y' z')/(x^2 + y^2))) +
BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x (x')^2)/
Sqrt[x^2 +
y^2] + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 y x' y')/
Sqrt[x^2 +
y^2] - (a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x (y')^2)/
Sqrt[x^2 +
y^2] + (4 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 x (z')^2)/
Sqrt[x^2 + y^2]) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y (x')^2)/
Sqrt[x^2 + y^2]) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 x x' y')/
Sqrt[x^2 +
y^2] + (a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y (y')^2)/
Sqrt[
x^2 + y^2] + (4 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 y (z')^2)/
Sqrt[x^2 + y^2])),
BesselJ[4,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a b^5 x y (-3 x^2 +
y^2) (x')^2)/(2 (x^2 + y^2)^2) + (a b^5 y^2 (-3 x^2 +
y^2) x' y')/(x^2 +
y^2)^2 - (a b^5 x y (x^2 +
5 y^2) (y')^2)/(2 (x^2 + y^2)^2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b^5 x^2 (x^2 -
3 y^2) (x')^2)/(2 (x^2 + y^2)^2) + (a b^5 x y (x^2 -
3 y^2) x' y')/(x^2 +
y^2)^2 + (a b^5 (x^4 + 3 x^2 y^2 -
2 y^4) (y')^2)/(2 (x^2 + y^2)^2))) +
BesselJ[3,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((2 a b^4 Sqrt[b^2 + c^2] y (-3 x^2 + y^2) t' x')/(x^2 +
y^2)^(3/2) + (3 a b^4 y (-3 x^2 y + y^3) (x')^2)/(x^2 +
y^2)^(5/2) + (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) t' y')/(x^2 + y^2)^(3/
2) - (6 a b^4 x y (-3 x^2 + y^2) x' y')/(x^2 + y^2)^(5/
2) + (3 a b^4 (-2 x^4 + 3 x^2 y^2 + y^4) (y')^2)/(x^2 +
y^2)^(5/2) - (2 a b^4 c y (-3 x^2 + y^2) x' z')/(x^2 +
y^2)^(3/
2) + (2 a b^4 x (4 c y^2 +
Sqrt[b^2 + c^2] (x^2 + y^2)) y' z')/(x^2 + y^2)^(3/2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((2 a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) t' x')/(x^2 +
y^2)^(3/2)) - (3 a b^4 x y (x^2 - 3 y^2) (x')^2)/(x^2 +
y^2)^(5/2) + (2 a b^4 Sqrt[
b^2 + c^2] y (-3 x^2 + y^2) t' y')/(x^2 + y^2)^(3/
2) + (6 a b^4 x^2 (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(5/
2) + (3 a b^4 x y (5 x^2 + y^2) (y')^2)/(x^2 + y^2)^(5/
2) + (2 a b^4 c x (x^2 - 3 y^2) x' z')/(x^2 + y^2)^(3/
2) - (2 a b^4 y (-2 c (x^2 - y^2) +
Sqrt[b^2 + c^2] (x^2 + y^2)) y' z')/(x^2 + y^2)^(3/2))) +
BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x t' x')/
Sqrt[x^2 +
y^2] - (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(x^2 + y^2)^(3/
2) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 y t' y')/
Sqrt[
x^2 + y^2] + (2 a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 x' y')/(x^2 + y^2)^(3/
2) + (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(x^2 + y^2)^(3/
2) - (2 a c (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x x' z')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2]) (c^4 +
3 c^3 Sqrt[b^2 + c^2] - b^2 (b^2 + 2 (b^2 + c^2)) +
c^2 (-2 b^2 + 3 (b^2 + c^2)) +
c (-4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) y y' z')/
Sqrt[x^2 +
y^2] - (4 a b^2 (c + Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 +
y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y t' x')/
Sqrt[x^2 + y^2]) + (a (c + Sqrt[b^2 + c^2])^2 (-b^2 +
2 c^2 + 2 c Sqrt[b^2 + c^2]) y^2 (x')^2)/(x^2 + y^2)^(3/
2) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 x t' y')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 + y^2)^(3/
2) + (a (c +
Sqrt[b^2 + c^2])^2 (2 (b^2 + c^2) x^2 + (3 b^2 +
c^2) y^2 + c^2 (2 x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (2 x^2 + y^2)) (y')^2)/(x^2 +
y^2)^(3/
2) + (2 a c (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y x' z')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2]) (b^4 + c^4 +
3 c^3 Sqrt[b^2 + c^2] + 2 b^2 (b^2 + c^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) +
c (4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) x y' z')/
Sqrt[x^2 +
y^2] - (4 a b^2 (c + Sqrt[b^2 + c^2])^2 y^2 (z')^2)/(x^2 +
y^2)^(3/2))) +
BesselJ[2,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] (-((a b x y ((-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + (2 b^2 +
c^2)^2 y^2 + c^4 (x^2 + y^2) +
4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (x')^2)/(2 (x^2 +
y^2)^2)) - (a b y^2 ((-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + (2 b^2 +
c^2)^2 y^2 + c^4 (x^2 + y^2) +
4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 + (a b x y ((2 b^2 + c^2)^2 x^2 + (5 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (y')^2)/(2 (x^2 + y^2)^2) - (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' z')/(x^2 +
y^2) + (8 a b^3 (c + Sqrt[b^2 + c^2]) y y' z')/(x^2 +
y^2) - (2 a b^3 (-3 c + Sqrt[b^2 + c^2]) (c +
Sqrt[b^2 + c^2]) x y (z')^2)/(x^2 + y^2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((a b x^2 (c^4 x^2 + (-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 +
y^2)) (x')^2)/(2 (x^2 +
y^2)^2)) - (a b x y (c^4 x^2 + (-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 -
1/(2 (x^2 + y^2)^2) a b ((2 b^2 + c^2)^2 x^4 + (3 b^4 +
2 b^2 (b^2 + c^2) + 3 (b^2 + c^2)^2) x^2 y^2 +
2 (-b^4 + (b^2 + c^2)^2) y^4 +
4 c Sqrt[
b^2 + c^2] (x^2 + y^2) ((2 b^2 + c^2) x^2 +
2 (b^2 + c^2) y^2) +
2 c^2 (x^2 + y^2) ((b^2 + 3 (b^2 + c^2)) x^2 +
6 (b^2 + c^2) y^2) + c^4 (x^4 + 3 x^2 y^2 + 2 y^4) +
4 c^3 Sqrt[
b^2 + c^2] (x^4 + 3 x^2 y^2 +
2 y^4)) (y')^2 + (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (x^2 - y^2) t' z')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) x y' z')/(x^2 +
y^2) - (2 a b^3 (c +
Sqrt[b^2 + c^2]) (-Sqrt[b^2 + c^2] x^2 +
c (x^2 - 2 y^2)) (z')^2)/(x^2 + y^2))) +
BesselJ[0,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(2 (x^2 +
y^2)) + (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y^2 x' y')/(x^2 +
y^2) - (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(2 (x^2 +
y^2)) + (2 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 + y^2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 (x')^2)/(2 (x^2 +
y^2)) + (a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 +
y^2) + (a b (c + Sqrt[b^2 + c^2])^2 ((3 b^2 + c^2) x^2 +
2 (b^2 + c^2) y^2 + c^2 (x^2 + 2 y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + 2 y^2)) (y')^2)/(2 (x^2 +
y^2)) +
4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' z' - (2 a b (c +
Sqrt[b^2 + c^2])^2 (b^2 x^2 + 2 c^2 (x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + y^2)) (z')^2)/(x^2 + y^2))),
BesselJ[4,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((a b^5 x y (5 x^2 +
y^2) (x')^2)/(2 (x^2 + y^2)^2) - (a b^5 x^2 (x^2 -
3 y^2) x' y')/(x^2 +
y^2)^2 - (a b^5 x y (x^2 -
3 y^2) (y')^2)/(2 (x^2 + y^2)^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((a b^5 (-2 x^4 + 3 x^2 y^2 +
y^4) (x')^2)/(2 (x^2 + y^2)^2)) - (a b^5 x y (-3 x^2 +
y^2) x' y')/(x^2 +
y^2)^2 - (a b^5 y^2 (-3 x^2 +
y^2) (y')^2)/(2 (x^2 + y^2)^2))) +
BesselJ[3,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((2 a b^4 Sqrt[b^2 + c^2] y (-3 x^2 + y^2) t' x')/(x^2 +
y^2)^(3/
2) - (3 a b^4 (x^4 + 3 x^2 y^2 - 2 y^4) (x')^2)/(x^2 +
y^2)^(5/2) + (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) t' y')/(x^2 + y^2)^(3/
2) + (6 a b^4 x y (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(5/
2) - (3 a b^4 x^2 (x^2 - 3 y^2) (y')^2)/(x^2 + y^2)^(5/
2) + (2 a b^4 y (4 c x^2 +
Sqrt[b^2 + c^2] (x^2 + y^2)) x' z')/(x^2 + y^2)^(3/
2) - (2 a b^4 c x (x^2 - 3 y^2) y' z')/(x^2 + y^2)^(3/
2)) + Cos[
Sqrt[b^2 + c^2] t -
c z] ((2 a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) t' x')/(x^2 +
y^2)^(3/2) - (3 a b^4 x y (x^2 + 5 y^2) (x')^2)/(x^2 +
y^2)^(5/2) - (2 a b^4 Sqrt[
b^2 + c^2] (-3 x^2 y + y^3) t' y')/(x^2 + y^2)^(3/
2) - (6 a b^4 y (-3 x^2 y + y^3) x' y')/(x^2 + y^2)^(5/
2) + (3 a b^4 x y (-3 x^2 + y^2) (y')^2)/(x^2 + y^2)^(5/
2) - (2 a b^4 x (2 c (x^2 - y^2) +
Sqrt[b^2 + c^2] (x^2 + y^2)) x' z')/(x^2 + y^2)^(3/
2) + (2 a b^4 c (-3 x^2 y + y^3) y' z')/(x^2 + y^2)^(3/
2))) + BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((2 a Sqrt[b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 y t' x')/
Sqrt[x^2 +
y^2] - (a (c + Sqrt[b^2 + c^2])^2 ((3 b^2 + c^2) x^2 +
2 (b^2 + c^2) y^2 + c^2 (x^2 + 2 y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + 2 y^2)) (x')^2)/(x^2 +
y^2)^(3/2) - (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x t' y')/
Sqrt[x^2 +
y^2] + (2 a (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 + y^2)^(3/
2) - (a (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 (y')^2)/(x^2 + y^2)^(3/

2) - (2 a (c + Sqrt[b^2 + c^2]) (b^4 + c^4 +
3 c^3 Sqrt[b^2 + c^2] + 2 b^2 (b^2 + c^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) +
c (4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) y x' z')/
Sqrt[x^2 +
y^2] + (2 a c (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y' z')/
Sqrt[x^2 +
y^2] + (4 a b^2 (c + Sqrt[b^2 + c^2])^2 x^2 (z')^2)/(x^2 +
y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((2 a Sqrt[b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 x t' x')/
Sqrt[x^2 +
y^2] - (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(x^2 + y^2)^(3/
2) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y t' y')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y^2 x' y')/(x^2 + y^2)^(3/
2) + (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(x^2 + y^2)^(3/
2) - (2 a (c + Sqrt[b^2 + c^2]) (c^4 +
3 c^3 Sqrt[b^2 + c^2] - b^2 (b^2 + 2 (b^2 + c^2)) +
c^2 (-2 b^2 + 3 (b^2 + c^2)) +
c (-4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) x x' z')/
Sqrt[x^2 +
y^2] - (2 a c (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y y' z')/
Sqrt[x^2 +
y^2] + (4 a b^2 (c + Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 +
y^2)^(3/2))) +
BesselJ[2,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] (-((a b x y ((5 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + (2 b^2 +
c^2)^2 y^2 + c^4 (x^2 + y^2) +
4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (x')^2)/(2 (x^2 +
y^2)^2)) + (a b x^2 ((2 b^2 + c^2)^2 x^2 + (-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 + (a b x y ((2 b^2 + c^2)^2 x^2 + (-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (y')^2)/(2 (x^2 + y^2)^2) + (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' z')/(x^2 +
y^2) + (8 a b^3 (c + Sqrt[b^2 + c^2]) x x' z')/(x^2 +
y^2) + (2 a b^3 (-3 c + Sqrt[b^2 + c^2]) (c +
Sqrt[b^2 + c^2]) x y (z')^2)/(x^2 + y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (1/(2 (x^2 +
y^2)^2) a b (-2 (b^4 - (b^2 + c^2)^2) x^4 + (3 b^4 +
2 b^2 (b^2 + c^2) +
3 (b^2 + c^2)^2) x^2 y^2 + (2 b^2 + c^2)^2 y^4 +
4 c Sqrt[
b^2 + c^2] (x^2 +
y^2) (2 (b^2 + c^2) x^2 + (2 b^2 + c^2) y^2) +
2 c^2 (x^2 +
y^2) (6 (b^2 + c^2) x^2 + (b^2 + 3 (b^2 + c^2)) y^2) +
c^4 (2 x^4 + 3 x^2 y^2 + y^4) +
4 c^3 Sqrt[
b^2 + c^2] (2 x^4 + 3 x^2 y^2 +
y^4)) (x')^2 + (a b x y ((-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + c^4 y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 + (a b y^2 ((-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + c^4 y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 +
y^2)) (y')^2)/(2 (x^2 + y^2)^2) - (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (-x^2 +
y^2) t' z')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) y x' z')/(x^2 +
y^2) + (2 a b^3 (c +
Sqrt[b^2 + c^2]) (-Sqrt[b^2 + c^2] y^2 +
c (-2 x^2 + y^2)) (z')^2)/(x^2 + y^2))) +
BesselJ[0,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(2 (x^2 +
y^2)) - (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 x' y')/(x^2 +
y^2) - (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(2 (x^2 +
y^2)) - (2 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 + y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((a b (c +
Sqrt[b^2 + c^2])^2 (2 (b^2 + c^2) x^2 + (3 b^2 +
c^2) y^2 + c^2 (2 x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (2 x^2 + y^2)) (x')^2)/(2 (x^2 +
y^2))) - (a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 +
y^2) - (a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y^2 (y')^2)/(2 (x^2 + y^2)) -
4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' z' + (2 a b (c +
Sqrt[b^2 + c^2])^2 (b^2 y^2 + 2 c^2 (x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + y^2)) (z')^2)/(x^2 + y^2))),
BesselJ[3,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a b^4 x (2 Sqrt[b^2 + c^2] (x^2 - y^2) +
c (x^2 + y^2)) (x')^2)/(x^2 + y^2)^(3/2) - (2 a b^4 Sqrt[
b^2 + c^2] y (-3 x^2 + y^2) x' y')/(x^2 + y^2)^(3/
2) + (a b^4 x (4 Sqrt[b^2 + c^2] y^2 +
c (x^2 + y^2)) (y')^2)/(x^2 + y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b^4 y (4 Sqrt[b^2 + c^2] x^2 +
c (x^2 + y^2)) (x')^2)/(x^2 + y^2)^(3/2) - (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(3/
2) + (a b^4 y (2 Sqrt[b^2 + c^2] (-x^2 + y^2) +
c (x^2 + y^2)) (y')^2)/(x^2 + y^2)^(3/2))) +
BesselJ[0,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] (4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' x' + (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x^2 x' z')/(x^2 +
y^2) + (4 a b^3 (c + Sqrt[b^2 + c^2])^2 x y y' z')/(x^2 +
y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' y' - (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y x' z')/(x^2 +
y^2) - (4 a b^3 (c + Sqrt[b^2 + c^2])^2 y^2 y' z')/(x^2 +
y^2))) +
BesselJ[2,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (x^2 - y^2) t' x')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) y (-x^2 +
y^2) (x')^2)/(x^2 + y^2)^2 + (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' y')/(x^2 +

y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) x (x^2 -
3 y^2) x' y')/(x^2 +
y^2)^2 - (16 a b^3 (c +
Sqrt[b^2 + c^2]) x^2 y (y')^2)/(x^2 +
y^2)^2 - (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x^2 x' z')/(x^2 +
y^2) - (4 a b^3 (c + Sqrt[b^2 + c^2])^2 x y y' z')/(x^2 +
y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' x')/(x^2 +
y^2)) - (16 a b^3 (c +
Sqrt[b^2 + c^2]) x y^2 (x')^2)/(x^2 +
y^2)^2 - (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (-x^2 +
y^2) t' y')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) y (-3 x^2 +
y^2) x' y')/(x^2 +
y^2)^2 + (8 a b^3 (c + Sqrt[b^2 + c^2]) x (-x^2 +
y^2) (y')^2)/(x^2 +
y^2)^2 + (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y x' z')/(x^2 +
y^2) + (4 a b^3 (c + Sqrt[b^2 + c^2])^2 y^2 y' z')/(x^2 +
y^2))) +
BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a (c + Sqrt[b^2 + c^2]) x (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (7 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
2 b^2 (-(b^2 - 2 (b^2 + c^2)) x^2 + (b^2 +
2 (b^2 + c^2)) y^2)) (x')^2)/(x^2 + y^2)^(3/
2) + (2 a (c + Sqrt[b^2 + c^2]) y (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (5 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
b^2 ((-3 b^2 + 2 (b^2 + c^2)) x^2 + (b^2 +
2 (b^2 + c^2)) y^2)) x' y')/(x^2 + y^2)^(3/
2) - (a (c + Sqrt[b^2 + c^2]) x (4 b^4 y^2 +
c^4 (x^2 + y^2) + 3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (3 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) (y')^2)/(x^2 +
y^2)^(3/2) + (8 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 x t' z')/
Sqrt[x^2 +
y^2] + (8 a b^2 (c + Sqrt[b^2 + c^2])^2 x y x' z')/(x^2 +
y^2)^(3/
2) - (8 a b^2 (c + Sqrt[b^2 + c^2])^2 x^2 y' z')/(x^2 +
y^2)^(3/2) - (4 a b^2 c (c + Sqrt[b^2 + c^2])^2 x (z')^2)/
Sqrt[x^2 + y^2]) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((a (c + Sqrt[b^2 + c^2]) y (4 b^4 x^2 +
c^4 (x^2 + y^2) + 3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (3 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) (x')^2)/(x^2 +
y^2)^(3/2)) + (2 a (c +
Sqrt[b^2 + c^2]) x (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (5 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
b^2 ((b^2 + 2 (b^2 + c^2)) x^2 + (-3 b^2 +
2 (b^2 + c^2)) y^2)) x' y')/(x^2 + y^2)^(3/
2) + (a (c + Sqrt[b^2 + c^2]) y (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (7 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
2 b^2 ((b^2 + 2 (b^2 + c^2)) x^2 - (b^2 -
2 (b^2 + c^2)) y^2)) (y')^2)/(x^2 + y^2)^(3/
2) + (8 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 y t' z')/
Sqrt[x^2 +
y^2] + (8 a b^2 (c + Sqrt[b^2 + c^2])^2 y^2 x' z')/(x^2 +
y^2)^(3/
2) - (8 a b^2 (c + Sqrt[b^2 + c^2])^2 x y y' z')/(x^2 +
y^2)^(3/2) - (4 a b^2 c (c + Sqrt[b^2 + c^2])^2 y (z')^2)/
Sqrt[x^2 + y^2]))};


To find the solution, I use the following procedure. First I take the values of the parameters to be:

b = 1
c = 20
a = 10^(-6)


Then I numerically evaluate them:

ngeodesic = N[geodesic];



And use the following procedure in order to obtain a solution for all the coordinates. First I calculate the first component of the four-velocity at parameter value equal to zero - this gives the remaining initial value for the NDSolve. After the NDSolve piece of code, I list the four-velocity vector's length to check if the system evolves the way it should.

minkowskimetric = {{-1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0,
0, 1}}

computeSoln[maxτi_, ivsi_, icsi_] :=
Block[{ivs, ics, i, χ, tmp, soln},
ics = icsi;
ivs = Join[{χ}, ivsi];
tmp = minkowskimetric;
tmp = ivs.(tmp.ivs);
χslv = Solve[tmp == uinvar, χ];
ivs[[1]] = Last[χ /. χslv];
deq = Table[
coords[[i]]''[τ] == (ngeodesic[[i]] /.
Join[Table[coords[[i]]' -> coords[[i]]'[τ], {i, 1, n}],
Table[coords[[i]] -> coords[[i]][τ], {i, 1, n}]]), {i, 1,
n}];
deq = Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}],
Table[coords[[i]][0] == ics[[i]], {i, 1, n}]];
soln = NDSolve[deq, coords, {τ, 0, maxτi}];
soln]
uinvar = -1;

udotu[solni_, τval_] := Block[{xα, uα},
xα =
Table[coords[[i]][τ] /. solni, {i, 1, n}] // Flatten;
uα = D[xα, τ];
xα = xα /. τ -> τval;
uα = uα /. τ -> τval;
uα.(minkowskimetric.uα)]
coordlist[τin_] :=
Table[ToString[coords[[i]]] <>
"=" <> {ToString[coords[[i]][τin] /. soln // First]}, {i, 1,
n}]


In what follows I use the initial values that actually yield some result. The behaviour looks like something obtained by my advisor, who used Maple. The ivs table is the initial velocities table - (x,y,z) components. The ics table is the table of initial coordinates (t,x,y,z).

maxτ = 40
ivs = {0.01, 0.01, 0.1};
ics = {0, 0.05, 0.05, 0.0};
soln = computeSoln[maxτ, ivs, ics];
ParametricPlot3D[
Evaluate[ {coords[[2]][τ], coords[[3]][τ],
coords[[4]][τ]}] /. soln, {τ, 0, maxτ}]
Plot[Evaluate[Table[coords[[i]][τ] /. soln, {i, 1, n}]], {τ,
0, maxτ}, PlotStyle -> Table[Hue[(i - 1)/n], {i, 1, n}],
PlotLegends -> coords, AxesLabel -> {"τ", "Coordinate"}]
ParametricPlot[
Evaluate[ {coords[[2]][τ], coords[[3]][τ]}] /.
soln, {τ, 0, maxτ}]

Join[{{"", "", "", "u*u values"}},
Table[{"τ=", ToString[i], "->", udotu[soln, i]}, {i, 0,
maxτ, maxτ/5}]] // TableForm


The problem begins when I try to change the initial conditions to, let's say:

ics = {0, 0.5, 0.05, 0.0};


In that case (and literally any other) I get the error:

NDSolve::ndsz: At τ == 0.34526573966566704, step size is effectively zero; singularity or stiff system suspected.


My advisor (who does not use Mathematica) gets very nice results for almost any imaginable parameters - since the physical problems only occur on the z coordinate axis. I am at loss at the moment - I suspect that the very small values (e.x. the amplitude "a") are the cause for the problem. However, I have no idea if I could even rescale all the quantities used here - the problem is, after all, highly nonlinear, even in the parameters a,b,c.

EDIT:

I put the whole code here: besselgravcart

besselgravcylindrical

besselgravconstants

• With default options, NDSolve becomes numerically unstable for ics[[2]] larger than about 0.1189. Method->"StiffnessSwitching" does not help. Neither does increasing WorkingPrecision (with all inputs rationalized). Perhaps, better choices of other options would. By the way, comments by other experienced users suggest that Maple is better than Mathematica at solving some ODEs. – bbgodfrey Feb 27 '19 at 23:41
• Does a symbolic expression exist for the inner turning point of the solution in {coords[[2]][τ], coords[[3]][τ]} space? Alternatively, does the solution satisfy any conservation laws? – bbgodfrey Feb 27 '19 at 23:55
• What do you mean by "a symbolic expression for the inner turning point of the solution"? I do not quite understand. Thanks for the question about the conservation laws! That might be an excellent way to help the stability. – K.T. Feb 28 '19 at 20:25
• Solutions in x - y space come in from large r, skim a circle of constant r for some distance and then go out to large r again. So, I speculate that this circle of constant r is a separatrix. If so, then infinitesimal errors in the computation near the separatrix would cause the computation to fail catastrophically. First integrals of the motion might help to alleviate this problem. Can the equations be rewritten conveniently in terms or r and theta instead of x and y (i.e., in cylindrical coordinates)? In other words, is the metric axisymmetric? – bbgodfrey Feb 28 '19 at 20:53
• Do not write part of your question as an answer. It will be deleted. Instead, post the lengthy code on the internet and provide a reference in the question. Also, your last comment tells me that t is a linear function of z, so one of the four ODEs is not needed. Correct? Is the metric also axisymmetric, which could reduce another equation to first order or eliminate it entirely? By the way, how many cycles in x - y do you need to compute to obtain a useful answer? – bbgodfrey Feb 28 '19 at 21:43