Stability-preserving positive neural ODE step

Integrate this into a positive neural ODE, recurrent state-space layer, or continuous-depth residual block whose hidden state is constrained to $u\geq0$. Let the network produce a vector field $f_\theta(u,x,t)$ and use the nonstandard update instead of Euler. First estimate or specify a damping bound $c_i>0$ for each state channel. For a softplus-parameterized positive state, estimate $c_i$ from a running maximum of $[-f_i(u)/\max(u_i,\epsilon)]_+$. Choose $\phi_i=h+\alpha_i h^2$ and compute $\Psi_i=\phi_i/(1+\tau_i\phi_i)$ with $\tau_i=\alpha_i$. Compute $f=f_\theta(u)$ and update $u_i\leftarrow u_i+\Psi_i f_i$. Before accepting the step, reduce $h$ until $\Psi_i\leq1/(c_i+\epsilon)$ for every channel. Every few training steps, estimate the largest eigenvalue of $D_hJ_f$ using three to five Jacobian-vector products and power iteration, then reduce $h$ or increase $\tau_i$ if $\rho(I+D_hJ_f)\geq1$. The paper supplies the denominator construction, positivity inequality, and Schur/Jury stability target; the neural implementation estimates $c_i$, Jacobian spectra, and local Taylor coefficients empirically. First test a four-to-eight-block positive neural ODE on MNIST or CIFAR-10, comparing equal-depth Euler and explicit midpoint at equal function evaluations. Measure negative-state fraction, exploding trajectories, gradient-norm variance, validation accuracy, and wall-clock loss descent. Success is zero sign violations and fewer unstable runs while matching or improving accuracy at the same function-evaluation budget.