Branch-Free Matrix-Function Methods for Second-Order Linear ODEs near Critical Damping and Resonance

Classical formulas for second-order constant-coefficient ODEs split the solution into distinct-root, repeated-root, and resonant cases; near the transition points, otherwise correct formulas may suffer cancellation or small-denominator loss. We derive a branch-free matrix-function representation for ay″+by′+cy=f(t), a≠0 using two entire functions of the scaled discriminant for homogeneous propagation and shifted ϕ-functions for polynomial–exponential forcing. The resulting expressions remain analytically defined at critical damping and forcing resonance. On bounded finite horizons, we establish a representation-level perturbation bound that remains nonsingular as the discriminant tends to zero and a conditional backward-error interpretation for the augmented matrix exponential; these statements remove representation-induced singularities but do not imply unconditional stability for arbitrary scales or times. Analytic coefficient sensitivities remain finite at critical damping. Reproducible 80-digit transition tests, randomized transition-focused ensembles, and timing benchmarks quantify the accuracy and cost of the method. We also clarify how the scalar kernel can be used for decoupled or method-of-steps subproblems arising in inertial neural-network models, while genuinely coupled nonlinear or delayed systems require a block/generalized extension.