EQ =
{R''[l] == 0.9897 R[l], R'[0] == 0, R[0] == 1,
Z'[l] == Sqrt[1 - R'[l]^2], Z[0] == 0};
NDSolve[EQ, {R, Z}, {l, 0, 2}];
{z[u_], r[u_]} = {Z[u], R[u]} /. First[%];
ParametricPlot[{z[l], r[l]}, {l, 0, 2},
PlotStyle -> {Red, Thick},
AspectRatio -> Automatic,
GridLines -> Automatic]
Table[{l, z[l], r[l], r'[l]}, {l, 0, 2, .2}] // TableForm
The output above first plot seems to succeed for portions of z[l] and r[l] for as long as they are real, however, it does not plot on real axes with r[l] and r'[l] even if they are real.
So, what is a numerical workaround to plot them? Taking the real part with zero imaginary part is not so elegant an option I feel.
(As remote connected background I refer to Zeeman/ R. Thom's views, but in mathematics, perhaps, there is no such catastrophe. We know what happens in a geometric singularity.)


ReorChop, but in earlier versions, that's the solution you find in these links. – Michael E2 Jul 10 '15 at 15:29