Why is there an error about division by 0?

I'm converting all of my old Mathematica 7.0 codes to my new Mathematica 13.2, and I'm struggling with some issues I didn't encounter with version 7.0. The following code is giving me a headache, since it worked perfectly with version 7.0 (this is a stripped down version of my full code, which is much more complex).

The following MWE is producing my issue:

axe = {0, 1}; (* The magnetic axis *)

coords[s_] = {x[s], z[s]};

r[s_] = Norm[coords[s]];

(* Note:  u is the radius of a spherical hole inside a magnetized sphere *)
dipole[s_, u_] =  Piecewise[{
{0, (0 <= r[s] < u)}, {1, (r[s] >= u)}
}] (3 (axe . coords[s]) coords[s]/r[s]^5 - axe/r[s]^3);

(* Note: 0.75 is the magnetic susceptibility of the sphere.  Could be any number *)
Coef1[s_, u_] = Piecewise[{
{1, (0 <= r[s] < u)},
{3 (3 - 0.75), (u <= r[s] < 1)},
{9 (1 - 0.75) + 2 0.75^2 (1 - u^3), (r[s] >= 1)}
}];

Coef2[s_, u_] = Piecewise[{
{0, (0 <= r[s] < u)},
{- 3 0.75 u^3, (u <= r[s] < 1)},
{0.75 (3 - 0.75) (1 - u^3), (r[s] >= 1)}
}];

orientation[s_, u_] = Normalize[Coef1[s, u] axe + Coef2[s, u] dipole[s, u]];

nlines = 50;

curves[u_, n_] := curves[u, n] = NDSolve[{
x'[s] == {1, 0} . orientation[s, u],
z'[s] == {0, 1} . orientation[s, u],
x[0] == 0.15 (2 n - 1 - nlines)/2,
z[0] == -3
}, {x, z}, {s, 0, 50},
Method -> BDF, PrecisionGoal -> 6, StoppingTest -> z[s] > 3]

graph[u_, n_] := ParametricPlot[Evaluate[{x[s], z[s]} /. curves[u, n]], {s, 0, 50}]

Manipulate[Show[Table[graph[u, n], {n, 1, nlines}],
PlotRange -> {{-1, 1}, {-1, 1}},
ImageSize -> 600], {{u, 0.5, "Radius"}, 0, 1, 0.01}
]


The code is numerically solving a differential equation to draw magnetic field lines around and inside a magnetized sphere with a hole inside. The code works and the output is nice (well at least in my full version), but Mathematica gives me lots of error messages about a division by 0. I'm unable to find what is wrong with the code (which was working without troubles in version 7.0).

So what's wrong in this code? Where is that division by 0 and how can I fix it? The only place where I could think of a division by 0 is in the Normalize command, for the definition of the function orientation[s, u] (the field line orientation should be an unit vector). But then, the Piecewise functions are supposed to prevent any problem here (??).

• Maybe I am wrong, but can you actually use NDSolve here? x'[s] =={1,0} . orientation[s, u] has the non-numerical quantity u in it and the boundary conditions also depend on the non-numeric quantity n. NDSolve is only called with one independent variable s. Maybe use something like ParametricNDSolve[{x'[s] == {1, 0} . orientation[s, u], z'[s] == {0, 1} . orientation[s, u], x[0] == 0.15 (2 n - 1 - nlines)/2, z[0] == -3, WhenEvent[z[s] > 3, "StopIntegration"]}, {x, z}, {s, 0, 50}, {u, n}]?
– ydd
Commented Jul 27, 2023 at 17:20
• @ydd, it doesn't work with this option. It's actually worst!
– Cham
Commented Jul 27, 2023 at 17:36
• Ah ok. I did get a plot using ParametricNDSolve but I have no idea if it's correct. You will probably have to run it on your local device to actually see the plots wolframcloud.com/obj/dtrimas/Published/… But it sounds like this is not what you're looking for
– ydd
Commented Jul 27, 2023 at 17:48
• Try Method -> {BDF, "DiscontinuityProcessing" -> None}. Getting a real fix might take real work. I happy enough with this (or rather I accept being happy cuz I don't have time for the alternative). Commented Jul 27, 2023 at 20:59
• @MichaelE2, wow! It works great. Thanks a lot!
– Cham
Commented Jul 27, 2023 at 21:24

Try Method -> {BDF, "DiscontinuityProcessing" -> None}. Getting a real fix might take real work. I happy enough with this (or rather I accept being happy cuz I don't have time for the alternative)
It might be a bug. Like I said, it would take some work to track that down, I suspect. If you have time, consider reporting it to WRI. NDSolve[] sets up events to handle Piecewise[]. Ideally it improves the solution. Since it's done automatically, I don't think it should result in division by zero. But you could see what WRI says.