# Surviving numerical cusp catastrophe

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.)

• Related (duplicates?): (17202), (34365), (75405) -- In V10.1, everything works fine without using Re or Chop, but in earlier versions, that's the solution you find in these links. – Michael E2 Jul 10 '15 at 15:29

Is this what you want to plot?

 ParametricPlot[Chop @ {z[l], r[l]}, {l, 0, 2},
PlotStyle -> {Red, Thick},
AspectRatio -> Automatic,
GridLines -> Automatic]


ParametricPlot[Chop @ {r[l], r'[l]}, {l, 0, 2},
PlotStyle -> {Red, Thick},
AspectRatio -> Automatic,
GridLines -> Automatic]


• When the imaginary component is very small making it numeric real, why cannot it be made really real by automatic chopping off the imaginary part inside the software? – Narasimham Jul 10 '15 at 18:35
• As noted in the comments to your question, in V10.1, the explicit application of Chop is not needed, so the problem you encountered seems to be fixed. – m_goldberg Jul 10 '15 at 21:36
• thanks. ( no future cusp catastrophes !) – Narasimham Jul 11 '15 at 13:12