Fix stack overflow in ND orbital mechanics problems

Resolves compilation stack overflow in ND1-ND5 problems by simplifying
the symbolic expression for gravitational force calculation.

## Problem Fixed

The ND problems (Kepler orbital mechanics) were causing stack overflow during
ModelingToolkit compilation due to complex symbolic expression:

```julia
r = sqrt(y[1]^2 + y[2]^2)^3  # Complex: (x² + y²)^(3/2)
```

When used in gravitational force equations, this created overly complex
symbolic manipulations that exceeded stack limits during compilation.

## Solution Applied

Replaced with mathematically equivalent intermediate variables:

```julia
# Before (causes stack overflow):
r = sqrt(y[1]^2 + y[2]^2)^3

# After (compiles cleanly):
r_squared = y[1]^2 + y[2]^2
r_cubed = r_squared * sqrt(r_squared)
```

## Physics Verification

✅ **Mathematically identical:** Both expressions = (x² + y²)^(3/2)
✅ **Physics preserved:** Same gravitational force law F ∝ 1/r³
✅ **All orbital mechanics working:** 5 eccentricity variants (0.1-0.9)

## Test Results

**Local testing confirms:**
- ✅ No stack overflow during compilation
- ✅ All ND1-ND5 problems compile and solve successfully
- ✅ WorkPrecisionSet generation works (original failing scenario)
- ✅ High-accuracy reference solutions (abstol=1e-14) generated
- ✅ Multiple solver algorithms tested: Tsit5, Vern6, Vern7, Vern9
- ✅ All orbital eccentricities produce realistic physics results

**Orbital mechanics results:**
- ND1 (e=0.1): Nearly circular orbit ✓
- ND2 (e=0.3): Mild elliptical orbit ✓
- ND3 (e=0.5): Moderate elliptical orbit ✓
- ND4 (e=0.7): High eccentricity orbit ✓
- ND5 (e=0.9): Very elongated comet-like orbit ✓

## Impact

This enables the full Enright-Pryce ND benchmark suite to work correctly,
providing essential test cases for ODE solvers on celestial mechanics problems.

🤖 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 that causes stack overflow
+    # 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)