Unexpected/unphysical behavior of ODE near singularity

Question

I'm trying to solve an ODE that becomes stiff or singular at a point. The ODE and initial conditions are

$$\begin{align*} \frac{d^2 z}{dt^2} &= f(z) f'(z) \bigg(1-\frac{3 \eta^2}{2} f(z) \bigg)\\ \frac{dz}{dt}\Bigg|_{t=t_0} &= \sqrt{1-\eta^2}\\ z\big|_{t=t_0}&=0. \end{align*}$$

where $f(z) = 1+\frac{z^2}{l^2} - \frac{\mu^d z^d}{l^{2d}}$. Here, $\mu$, $d$, and $l$ are constants, $\eta$ is an adjustable parameter between $0$ and $1$, and $t_0$ denotes the initial time.

I've implemented the system in the code snippet below:

d = 3;
l = 5.; μ = 26.^(1/d);

f[z_] := 1 + z^2/l^2 - μ^d z^d/l^(2 d);

Subscript[θ, 0] = 2.4;
tSoln = ParametricNDSolve[{ZZ''[t] == f[ZZ[t]] f'[ZZ[t]] (1 - (3 η^2)/2 f[ZZ[t]]), 
    ZZ'[-l Subscript[θ, 0]] == Sqrt[1 - η^2], 
    ZZ[-l Subscript[θ, 0]] == 0}, {ZZ}, 
    {t, -l Subscript[θ, 0], l Subscript[θ, 0]}, {η}]
z[t_, η_] := ZZ[η][t] /. tSoln;

Plot[z[t, 0], {t, -5*2.4, 5*2.4}]

From the underlying physics, I know that for sufficiently small $\eta$, the function $z(t)$ should asymptotically approach a constant value (here $z \to 25$). However, I believe the system becomes stiff/singular as $z(t)$ approaches that value, and Mathematica doesn't find the correct solution.

As an example, the above code gives me this plot of $z(t)$:

For some reason, the graph dips downward. (With other small, nonzero values of $\eta$, the graph crosses the asymptote of $z=25$.) If I take the $\eta=0$ case and integrate it directly rather than via an ODE solver, I obtain the behavior that I expect:

Is there an ODE solution method that is particularly suited for this problem? I've mixed and matched with "BDF", "StiffnessSwitching", and some others without really understanding how they work, and none of them give consistent results. How should I approach this?

In short: How can I get Mathematica to understand the correct asymptotic behavior of a differential equation when it becomes stiff or singular?

Sorry for such a lengthy post, it's kind of hard for me to articulate the exact problem. I'm also not very familiar with the advanced tools of Mathematica, so I very much appreciate any guidance. Thank you so much! (Please let me know if I can clarify anything.)

Answer 1

It's just a matter of precision. WorkingPrecision -> 32 will resolve your problem:

d = 3;
l = 5; μ = 26^(1/d);

f[z_] = 1 + z^2/l^2 - μ^d z^d/l^(2 d);

Subscript[θ, 0] = 24/10;
tSoln = ParametricNDSolve[{ZZ''[t] == f[ZZ[t]] f'[ZZ[t]] (1 - (3 η^2)/2 f[ZZ[t]]), 
    ZZ'[-l Subscript[θ, 0]] == Sqrt[1 - η^2], 
    ZZ[-l Subscript[θ, 0]] == 0}, {ZZ}, {t, -l Subscript[θ, 0], 
    l Subscript[θ, 0]}, {η}, WorkingPrecision -> 32];
z[t_, η_] = ZZ[η][t] /. tSoln;

Plot[z[t, 0], {t, -5*2.4, 5*2.4}]