Fix stack overflow in ND problems by simplifying symbolic expressions

The ND problems (ND1-ND5) in the Enright-Pryce benchmark suite were causing
stack overflow errors during ModelingToolkit compilation. The issue was caused
by the complex symbolic expression `r = sqrt(y[1]^2 + y[2]^2)^3` which created
overly complex symbolic manipulations when used in the gravitational force equations.

This fix breaks down the calculation into intermediate variables:
- r_squared = y[1]^2 + y[2]^2
- r_cubed = r_squared * sqrt(r_squared)

This mathematically equivalent formulation is much simpler for the symbolic
system to handle, resolving the compilation stack overflow while maintaining
the same physics (Kepler orbital dynamics).

Fixes the stack overflow issues reported in the ND benchmark suite.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

benchmarks/NonStiffODE/enright_pryce.jl

@@ -502,11 +502,14 @@ const NC_PROBLEMS = [nc1prob, nc2prob, nc3prob, nc4prob, nc5prob]
 y = collect(y)
 @parameters ε
 nd1sys = complete(let
-    r = sqrt(y[1]^2 + y[2]^2)^3
+    # Use intermediate variable to avoid complex symbolic expansion
+    # r_cubed = (x^2 + y^2)^(3/2) for the gravitational force law
+    r_squared = y[1]^2 + y[2]^2
+    r_cubed = r_squared * sqrt(r_squared)
     nd1eqs = [D(y[1]) ~ y[3],
         D(y[2]) ~ y[4],
-        D(y[3]) ~ (-y[1]) / r,
-        D(y[4]) ~ (-y[2]) / r
+        D(y[3]) ~ (-y[1]) / r_cubed,
+        D(y[4]) ~ (-y[2]) / r_cubed
     ]
     ODESystem(nd1eqs, t, name = :nd1)
 end)