Add higher-order differential equation support

[KingArth0r]

Summary

• add symbolic higher-order homogeneous constant-coefficient DSolve support
• add numeric higher-order IVP support by reducing explicit ODEs to RK4 systems
• cover higher-order symbolic and numeric differential equations in tests

Validation

• npm run lint
• node_modules\.bin\jest.cmd --config ./config/jest.config.cjs test/compute-engine/differential-equations.test.ts --runInBand --reporters default
• npx tsc --target es2022 --module es2022 --moduleResolution bundler --types node --skipLibCheck --allowImportingTsExtensions true --noEmit --ignoreConfig ./src/compute-engine.ts
• npm run check:deps

[arnog]

Excellent PR. Looks good overall, ready to merge. A few points for consideration and potential improvement are noted below.

1/ Integration-constant naming is inconsistent between orders. First-order still emits single capitals
via integrationConstant (C, then K, L, …), while second/higher-order emit c_1, c_2, … via the new
integrationConstants. So y'=y → C·e^x but y''=y → c_1·e^x + c_2·e^{-x}. Within one engine the user sees
two conventions. Recommend unifying on one (e.g. route first-order through c_1 too, or continue the
capital sequence for higher order). Minor: toString renders the new symbols quoted ("c_1"); worth
confirming they serialize to LaTeX as a clean c_1 subscript (they should).

2/ integrationConstant candidate list contains reserved/operator names. The new list includes D, E, N.
E (Euler) is filtered because it's a value-def with a declared type, but D (derivative) and N (numeric
eval) are operator defs, so isValueDef is false and they could be selected as a constant symbol in edge
cases (only reached once C,K,L,M,N,A,B are all taken; rare, but it'd produce confusing output like y =
D·e^x). Also the binding lookup only checks ce.context.lexicalScope.bindings (current scope), so
globally-bound names may not be filtered. Consider excluding any name that resolves to a binding
(operator or value).

3/ Order-≥3 exactness silently degrades to floats on complex roots. polynomialRoots returns only real
roots at degree ≥3 (complex ones are dropped), so whenever the characteristic polynomial has a complex
root the total-multiplicity check fails and the code falls to the numeric
solveHigherOrderWithNumericRoots path, emitting float coefficients. Visible in the PR's own test: y'''=y → cos(0.8660254…·x) instead of cos(\tfrac{\sqrt3}{2}x). It only looks exact when numeric roots land on
clean values (that's why y''''−y=0 shows exact cos x/sin x). Not wrong, just not exact where an exact
form exists. Worth a note in the code/PR that ≥3 with complex roots is numeric, and a future improvement
would be to recover exact complex roots (or chop-recognize √3/2).

4/ Maintainability: two near-parallel solver paths. solveSecondOrderHomogeneousConstantCoefficient
(~70 lines) and solveHigherOrderHomogeneousConstantCoefficient (~90 lines) overlap substantially. The
order-2 special case exists to get exact radical/complex forms via the quadratic formula (a real benefit,
given (3)), so keeping it is defensible but a one-line comment explaining why order 2 doesn't route
through the general path would help for future maintenance.

5/ Latent (pre-existing) derivative-form gaps, natural to close here. Order detection
(derivativeOrderOfDependent) handles nested D(D(y,x),x) and Apply(Derivative(f,n),x), but not the flat
multi-variable form D(y,x,x), which CE also produces and preserves (to be fair I've just added it coincidentally today)
it's mis-counted as order 1, so a flat-form 2nd-order ODE would be silently mis-solved as first order. This predates the PR (the old
isDerivativeOfDependent had the same blind spot) and isn't reachable from y''(x) LaTeX (which parses to
the nested form), but the new order-counting code is the natural place to handle it (count the
independent-variable args). Separately, Leibniz \frac{d^2y}{dx^2} doesn't parse as a derivative at all; a pre-existing parser gap, out of scope here but worth logging.
a pre-existing parser gap, out of scope here but worth logging.

6/ rk4System finiteness handling. It early-returns undefined after each stage evaluation, whereas
scalar rk4 computes all four k's then checks once. The new behavior is arguably better (avoids
propagating NaN), just slightly inconsistent with its sibling, fine to leave, or align for symmetry.